TRCA-OSAP Stream Habitat Analysis

Introduction

Background

The ability of streams and rivers to support a healthy and diverse biota depends on the availability of suitable habitat. This includes channel structure (e.g., cover, substrate) and channel processes (e.g., hydrology and sediment transport). Understanding these habitat components of a river system is crucial as it influences habitat availability, water quality, and overall ecosystem health. In the context of Ontario Stream Assessment Protocol (OSAP) and the Toronto and Region Conservation Authority (TRCA), in-stream habitat monitoring can provide valuable insights into the health of aquatic ecosystems, particularly for fish and benthic invertebrates, which are often sensitive to changes in their physical environment. By examining the relationships between habitat conditions and biological responses, we can better understand the impacts of natural and anthropogenic disturbances on aquatic ecosystems to guide land-use planning and restoration efforts.

Project Overview

This study aims to understanding the in-stream habitat data collected through TRCA’s long-term monitoring programs. It includes compilation and summarization of the data, as well as exploring spatial and temporal patterns. The results will allow TRCA and OSAP to assess the variability in habitat across the region and how they relate to fish and benthic invertebrate communities to better inform monitoring and management decisions.

TRCA Monitoring Locations

Methods

Data and models

Habitat and fish community data were downloaded from the Flowing Waters Information System (FWIS) Channel Structure Summary table. Watershed outlines were downloaded from TRCA open data portal. Supplemental data was joined from V3 of the Ontario Aquatic Ecosystem Classification (AEC). The Upstream catchment stress index was used to generate 3 Stress Classes corresponding to low (20-80), medium (80-140), and high (140-200) anthropogenic impacts. The Upstream Catchment Area (UCA) was used to generate 3 UCA Classes corresponding to small (0-20), intermediate (20-80), and large (>80) stream sizes. Finally, the site median Habitat Stability Class from the OSAP data was used to generate 2 Habitat Stability Classes (> and <= 1) corresponding to low and high stability.

Generalized Linear Mixed Effects Models (GLMMs; using the GLMMtmb R package) were used to describe spatial and temporal patterns in in-stream habitat endpoints over time. Distributional families were selected to match the distributions and data generating processes that underlie the endpoints, for example substrate sizes tend to vary on a log-normal distribution, whereas percentages and cover endpoints generally follow a beta distribution.

Some endpoints contained multiple categorical values within continuously measured values (i.e., sand, silt, clay categories for fines within the Substrate Sizes endpoints). These values were all replaced by 0’s and modeled with the zero-inflation portion of the model. See the R Code for full list of data transformations. Endpoints with more than 75% of observations at 0 tended to have worse preliminary model fits and were excluded from the final analysis.

Each endpoint was modeled with the following formula: \[ Endpoint = Developed + Agricultural + ln(Reach Slope+1) + Year_n*UCA*HabitatStability + Year_n*UCA*Stress + (Year_n|SiteCode) \]

where \(Developed\) and \(Agricultural\) are percent in the upstream catchment (from AEC); \(Reach Slope\) is the percent slope of the reach from the AEC; \(UCA\), \(HabitatStability\), and \(Stress\) are fixed effect factor terms; \(Year_n\) represent temporal slopes and the \((Year_n|SiteCode)\) is a nested random-effects term which fits a temporal slope to each monitoring location nested within a UCA class, Habitat Stability class and Stress Class factor. The \(Year_n*UCA*HabitatStability\) and \(Year_n*UCA*Stress\) are both fit as three-way interaction terms fitting temporal slopes to the three factor variables. Data coverage was too small to fit a four-way interaction term (i.e., there were insufficient sites corresponding to each factor level). Some endpoints are also fit with a zero-inflation model (probability of a zero: \(E_0\)). Zero-inflation was modeled the same terms as above.

The following table summarizes the distributional family used for each endpoint.

	Distribution
Cover
Wood Cover - Unembedded	beta
Wood Cover - Total	beta
Round Rock Cover - Unembedded	beta
Round Rock Cover - Total	beta
Flat Rock Cover - Unembedded	beta
Flat Rock Cover - Total	beta
Percent
Filamentous Algae Percent	beta
Non-Filamentous Algae Percent	beta
Moss Percent	beta
Indices
Sediment Sorting Index	gamma
Bank Stability - Mean	beta
Habitat Homogeneity	beta
Substrate Sizes
D50	lognormal
D85	lognormal
D50Max	lognormal
Morphology
Width - Mean	lognormal
Depth - Mean	lognormal
Width:Depth Ratio	lognormal
Usable Area	gamma

After an initial round of model fitting, models with convergence errors were subject to an outlier removal method. 10,000 simulated residuals were generated from the fitted models, and observed residuals were compared to simulated residuals using the DHARMa R package. Observed residuals that were outside the 2.5th and 97.5th percentiles of simulated residuals were considered outliers. These values were removed from the model training data, and models were refit to the reduced dataset. Outlier values were excluded from model fitting, but included in final evaluations.

R Code

This section contains all the R code to reproduce the current results and analysis.

0 FWIS Source

if( !require( "httr" ) ) {
  install.packages( "httr" )
}
library( httr )
if( !require( "jsonlite" ) ) {
  install.packages( "jsonlite" )
}
library( jsonlite )
if( !require( "base64enc" ) ) {
  install.packages( "base64enc" )
}
library( base64enc )

# widePost.r -- a collection of functions that enable access to WIDE server data using R
#
# Usage: source( "widePost.r" )  ...from R console/script, with current directory set to folder with widePost.r
#
# to enable some debugging, such as the HTTP request/response headers and JSON request, enter the following:
#  options( verbose = 1 )
# to disable subsequent debugging, use:
#  options( verbose = 0 )

.WIDE <- new.env()
.WIDE$clientVersionDate <- "2016-06-27-1400"
.WIDE$clientVersion <- c( "1", "0" )
# list of servers that the user has logged in to...
.WIDE$connectedServers <- list()

# wideLogin -- connect to a WIDE data API server
#  serverURL: the URL of the server data API script
#  username/password: login credentials
# returns false in the event of an error, or a serverHandle (integer > 0) value to be passed to subsequent widePost functions

wideLogin<-function( serverURL, username, password )
{
  # application.php requires parameters "Username", "Password" and "cls_logon_submit" as follows...
  postBody<-list( Username = username, Password = password, cls_logon_submit = 1 )
  loginResponse <- httr::POST( serverURL, body = postBody, encode = "form" )
  loginContent <- content( loginResponse, "text" )
  
  # error handling
  if( !jsonlite::validate( loginContent ) ) {
    # php syntax errors do not get encoded as JSON, so check for beginning and ending curly braces...
    if( ( substr( loginContent, 1, 1 ) != "{" ) &&
        ( substr( loginContent, nchar( loginContent ), nchar( loginContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", loginContent ) )
    } else {
      print( paste( " invalid JSON: ", loginContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( loginContent )
  
  if( "error" %in% names( jsonData ) ) {
    print( jsonData )
    return( FALSE )
  }
  
  # successful login
  print( jsonData$status )
  
  # check for version differences
  if( .WIDE$clientVersion[ 1 ] == jsonData$version[ 1 ] ) {
    if( .WIDE$clientVersion[ 2 ] == jsonData$version[ 2 ] ) {
      print( paste( "(0) WIDE interface ready - client/server version", paste( .WIDE$clientVersion, collapse = "." ),
                    " (", .WIDE$clientVersionDate, ")" ) )
      versionStatus <- 0
    } else if( .WIDE$clientVersion[ 2 ] > jsonData$version[ 2 ] ) {
      versionStatus <- 1
      print( paste0( "(1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                     " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (",
                     .WIDE$clientVersionDate, ")" ) )
    } else {
      versionStatus <- -1
      cat( paste0( "(-1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (",
                   .WIDE$clientVersionDate, ")\nPlease consider downloading a new client\n" ) )
    }
  } else if( .WIDE$clientVersion[ 1 ] > jsonData$version[ 1 ] ) {
    versionStatus <- 1
    print( paste0( "(1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (", .WIDE$clientVersionDate, ")" ) )
  } else {
    versionStatus <- -1
    print( paste0( "(-1) FATAL ERROR: client is out of date - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (", .WIDE$clientVersionDate,
                   ")\nPlease consider downloading a new client\n" ) )
  }
  
  # create and store server handle to be returned to the user
  if( versionStatus == -1 ) {
    return( -1 )
  }
  .WIDE$connectedServers <- append( .WIDE$connectedServers, serverURL )
  return( length( .WIDE$connectedServers ) )
}

# wideDataSelect -- retrieve one or more records from a database table
#  e.g.,
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName ASC"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName DESC"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list('"1" DESC'))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list('1 DESC'))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list(1))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),list(StreamName=c("like","%Black%")))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),list(StreamName=c("like","%Black%"),StreamCode=c("like","bl%")),list(1),"OR")
# returns a data.frame if successful; otherwise, returns NULL

wideDataSelect<-function( serverHandle, table, colsSelected, where = NULL, orderBy = NULL, logical = "AND", geoConstraint = 0 )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataSelect", serverHandle = serverHandle, table = table, colsSelected = colsSelected,
                 where = where, orderBy = orderBy, logical = logical, geoConstraint = geoConstraint ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: [", seq_len( length( .WIDE$connectedServers ) ), "]" ) )
  }
  
  # Note that the "xxx=" prefix in each "element" of the jsonPost list must be specified (i.e., it can't be left out
  #  of the jsonPost list() below)...
  # e.g., jsonPost<-list(table="lkpStreamCode",action="select",colsSelected=list("StreamName"))
  jsonPost<-list( table = table, action = "select", colsSelected = colsSelected, colsWhere = where, orderbyCols = orderBy, logical = logical, geoConstraint = geoConstraint )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  # for some circumstances (e.g., application development, testing, etc.) it may be useful to save the basic postResponse content
  #  to a file...
  # write( paste( "Content of postResponse ", content( postResponse, "text" ) ), "c:/temp/routput.txt" )
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( NULL )
  }
  
  jsonData <- fromJSON( postContent )
  
  if( "error" %in% names( jsonData ) ) {
    print( jsonData )
    return( NULL )
  }
  
  # find all of the binary columns which must be decoded
  binaryCols <- Filter( function( colName ) { jsonData[[ 1, colName ]] == 1 }, names( jsonData ) )
  
  # removing a row from a dataframe with one column ends up converting that data frame into a vector;
  # to prevent this from happening, use drop=FALSE
  jsonData <- jsonData[ -1, , drop = FALSE ]
  
  # if the data frame isn't empty after removing the row, renumber all of the rows
  if( nrow( jsonData ) ) {
    row.names( jsonData ) <- 1:nrow( jsonData )
  }
  
  # if there are any binary columns to decode, do that now... (base64decode cannot NA, so those must be explicitly
  #  returned back as NA)
  # NA_character_ is returned because that is the type of NA value used in the shinyapp DataTable
  if( length( binaryCols ) ) {
    jsonData[[ binaryCols ]] <- sapply( jsonData[[ binaryCols ]], function( i ) {
      if( is.na( i ) ) {
        return( NA_character_ )
      }
      # for some reason, without the paste (... collapse), the as.character(...) ends up as type list
      #  in the final jsonData dataframe
      return( paste0( "0x", as.character( base64decode( i ) ), collapse = "" ) )
    } )
  }
  
  return( jsonData )
}

# wideDataSpatialConstraint -- define a geospatial constraint for subsequent queries with
#   spatial constraints enabled
#  e.g.,

# return TRUE if successful; otherwise, return FALSE
wideDataSpatialConstraint <- function( serverHandle, table = NULL, column = NULL, where = NULL, comparisonResult = NULL, logical="AND") {
  
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataSpatialConstraint", serverHandle = serverHandle, table = table, column = column, where = where,
                 comparisonResult = comparisonResult, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "geospatialConstraint", column = column, colsWhere = where,
                  comparisonResult = comparisonResult, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# wideDataUpdate -- update one or more records of a database table
#  e.g.,
# wideDataUpdate(serverHandle = 1L, table = "lkpStreamCode", cols = list(StreamName = "Teststream1"), where =list(StreamCode = c("=", "teststr")))
# return TRUE if successful; otherwise, return FALSE

wideDataUpdate<-function( serverHandle, table, cols, where, logical = "AND" )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataUpdate", serverHandle = serverHandle, table = table, cols = cols, where = where, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  # Note that the "xxx=" prefix in each "element" of the jsonPost list must be specified (i.e., it can't be left out of the jsonPost list() below)...
  jsonPost<-list( table = table, action = "update", cols = cols, colsWhere = where, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  
  return( TRUE )
}

# return TRUE if successful; otherwise, return FALSE

wideDataInsert<-function( serverHandle, table, cols )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataInsert", serverHandle = serverHandle, table = table, cols = cols ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "insert", cols = cols )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid JSON: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# return TRUE if successful; otherwise, return FALSE

wideDataDelete<-function( serverHandle, table, where, logical = "AND" )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataDelete", serverHandle = serverHandle, table = table, where = where, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "delete", colsWhere = where, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid JSON: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# example:
#   data <- wideDataSelect("lkpStreamCode", list("*"), orderBy=list("StreamCode DESC"))
#   exportCSV(data, "c:/temp/output.csv")

exportCSV<-function( data, filepath )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "exportCSV", data = data, filepath = filepath ) )
  }
  
  # PARAMETERS (row.names, na)
  # - do not include row numbers
  # - save null values as conventional nulls instead of empty strings (i.e., "column1","column2",,"column4" instead of "column1","column2","","column4")
  write.csv( data, filepath, row.names = FALSE, na = "" )
}

# to read a csv, see read.csv(...)

1 Source

library(tidyverse)
library(sf)

# Columns that uniquely identify observations
ID_cols <- c(
  "ID",
  "StreamName",
  "StreamCode",
  "SiteCode",
  "Latitude",
  "Longitude",
  "SampleDate",
  "Sample",
  'OrgNameList'
)

# Endpoint names and groups
met_cols <- list(
  `Morphology` = list(
    `Width - Mean` = "AvgWettedWidth",
    `Depth - Mean` = "meanDepth",
    `Width:Depth Ratio` = "widthDepthRatio",
    `Usable Area` = "UsableArea"
  ),
  `Substrate Sizes` = list(
    `D15` = "Substd15",
    `D50` = "Substd50",
    `D85` = "Substd85",
    `D50Max` = "substPartMaxd50"
  ),
  `Indices` = list(
    #`Substrate Quality Rating` = "Substrate Quality Rating",
    `Sediment Sorting Index` = "sedimentSortingIndex",
    `Sediment Transport Potential` = "sedimentTransportPotential",
    `Bank Stability - Mean` = "bankStabilityAVG",
    `Habitat Homogeneity` = "scoreHomogeneity",
    `Habitat Stability Score` = "HabitatStabilityScore"
  ),
  `Percent` = list(
    `Filamentous Algae Percent` = "TotalVegTypePCFL",
    `Non-Filamentous Algae Percent` = "TotalVegTypePCAL",
    `Moss Percent` = "TotalVegTypePCSS",
    `Macrophytes Percent` = "TotalVegTypePCMC"
  ),
  `Cover` = list(
    `Wood Cover - Unembedded` = "UnembCoverTypePCWood",
    `Wood Cover - Total` = "TotalCoverTypePCWood",
    `Round Rock Cover - Unembedded` = "UnembCoverTypePCRndRock",
    `Round Rock Cover - Total` = "TotalCoverTypePCRndRock",
    `Flat Rock Cover - Unembedded` = "UnembCoverTypePCFltRock",
    `Flat Rock Cover - Total` = "TotalCoverTypePCFltRock",
    `Macrophyte Cover - Unembedded` = "UnembCoverTypePCMacrophyte",
    `Macrophyte Cover - Total` = "TotalCoverTypePCMacrophyte",
    `Bank Cover - Unembedded` = "UnembCoverTypePCBank",
    `Bank Cover - Total` = "TotalCoverTypePCBank"
  )
)

# Selected Fish Endpoints
sel_fish_ep<-c(
  `Community Total Biomass`="Comm_Biomass",
  `Community Total Density`="Comm_Abundance",
  `Community Total Richness`="Comm_Richness",
  `SATI Biomass`="Therm_SATI_Bioperc",
  `Cold Water Stenotherm Biomass`="Therm_cold_Bio",
  `% Cold Water Stenotherm Biomass`="Therm_cold_Bioperc",
  `Tolerant Taxa Biomass`="Tol_tol_Bio",
  `% Tolerant Taxa Biomass`="Tol_tol_Bioperc",
  `Intolerant Taxa Biomass`="Tol_intol_Bio",
  `% Intolerant Taxa Biomass`="Tol_intol_Bioperc",
  
  `Tolerant Taxa Richness`="Tol_tol_Rich",
  `Intolerant Taxa Richness`="Tol_intol_Rich",
  
  `Sensitive Darters Biomass`="Spec_DRT_Bio",
  `% Sensitive Darters Biomass`="Spec_DRT_Bioperc",
  `Brook Trout Biomass`="Spec_BKT_Bio",
  `% Brook Trout Biomass`="Spec_BKT_Bioperc",
  `Creek Chub Biomass`="Spec_CKC_Bio",
  `% Creek Chub Biomass`="Spec_CKC_Bioperc",
  `Open Substratum/Substratum Choosers Biomass`="Rep_OSS_Bio",
  `% Open Substratum/Substratum Choosers Biomass`="Rep_OSS_Bioperc",
  `Lithophilic Spawners Biomass`="Rep_Lith_Bio",
  `% Lithophilic Spawners Biomass`="Rep_Lith_Bioperc"
  
)

# function to reformat model terms

get_terms <- function(mod,add_fixed_var=T){
  fe_tbl <- as.data.frame(glmmTMB::fixef(mod)$cond)
  colnames(fe_tbl)<-"coef"
  fe_tbl$se <- diag(glmmTMB::sandwich(mod))
  
  add_var <- fe_tbl["(Intercept)","se"]
  if (!add_fixed_var) add_var <- 0
  
  re_tbl <- as.data.frame(glmmTMB::ranef(mod, condVar = T)) %>% 
    filter(component == "cond") %>% 
    mutate(val=condval,
           sd=condsd)
  
  comp_sel <- unique(re_tbl$grpvar)
  u_comp_sel <- str_split(comp_sel,"\\:")
  u_comp_sel <- unique(unlist(u_comp_sel))
  u_comp_sel <- unique(gsub("`","",u_comp_sel))
  
  for (i in u_comp_sel) {
    comp_sel_sub <- comp_sel[sapply(comp_sel,function(x) any(grepl(i,x)))]
    re_tbl[[i]]<-NA_character_
    
    for (ii in comp_sel_sub) {
      pos <- str_split(ii,"\\:",simplify = T)
      pos <- gsub("`","",pos)
      pos <- which(pos==i)
      
      re_tbl[re_tbl$grpvar==ii,i] <- str_split(re_tbl[re_tbl$grpvar==ii,"grp"],"\\:",simplify= T)[,pos]
    }
  }
  
  for (i in unique(re_tbl$SiteCode[!is.na(re_tbl$SiteCode)])){
    
    for (ii in unique(re_tbl$term)){
      val_sum <- re_tbl$val[re_tbl$SiteCode %in% i & re_tbl$term %in% ii]
      sd_sum <- re_tbl$sd[re_tbl$SiteCode %in% i & re_tbl$term %in% ii]
      
      for (iii in unique(re_tbl$grpvar)[!grepl("SiteCode",unique(re_tbl$grpvar))]){
        val_sum <- c(val_sum, re_tbl$condval[re_tbl$grpvar %in% iii &
                                               re_tbl$grp %in% gsub(paste0(":",i),"",unique(re_tbl[["grp"]][re_tbl$SiteCode %in% i])) &
                                               re_tbl$term %in% ii
        ]
        )
        sd_sum <- c(sd_sum, re_tbl$condsd[re_tbl$grpvar %in% iii &
                                            re_tbl$grp %in% gsub(paste0(":",i),"",unique(re_tbl[["grp"]][re_tbl$SiteCode %in% i])) &
                                            re_tbl$term %in% ii
        ]
        )
      }
      
      re_tbl$val[re_tbl$SiteCode %in% i & re_tbl$term %in% ii] <- sum(val_sum)
      re_tbl$sd[re_tbl$SiteCode %in% i & re_tbl$term %in% ii] <- sqrt(sum(sd_sum^2)+add_var^2)
      
    }
  }
  
  comp_list <- re_tbl %>% 
    filter(!is.na(SiteCode))%>% 
    mutate(group_var="SiteCode") %>% 
    select(group_var,group=SiteCode,component,term,val,sd,any_of(u_comp_sel)) %>% 
    group_by(group_var,group,component,term) %>% 
    summarise(
      value = sum(val),
      sd = sqrt(sum(sd^2)+add_var^2),
      across(any_of(u_comp_sel),~unique(.x)),
      .groups = "drop"
    )
  
  comp_list <- bind_rows(
    comp_list,
    re_tbl %>% 
      filter(is.na(SiteCode)) %>% 
      filter(!grpvar %in% unique(comp_list$group_var)) %>% 
      select(group_var=grpvar,
             group=grp,
             component,
             term,
             val=condval,
             sd=condsd,
             any_of(u_comp_sel)) %>% 
      group_by(across(c(any_of(u_comp_sel),group_var,group,component,term))) %>% 
      summarise(
        value = sum(val),
        sd = sqrt(sum(sd^2)+add_var^2),
        .groups = "drop"
      )
  )
  
  comp_list <- comp_list %>% 
    mutate(
      pred_low = value - (qnorm(1 - 0.1 / 2) * sd),
      pred_high = value + (qnorm(1 - 0.1 / 2) * sd),
      pred_low_trend = value - (qnorm(1 - 0.5 / 2) * sd),
      pred_high_trend = value + (qnorm(1 - 0.5 / 2) * sd),
    ) 
  
  return(comp_list)
  
}

get_terms2 <- function(mod,add_fixed_var=F,terms){

  fixed_list <- list(
    "`UCA Class`:`Stress Class`:`Habitat Stability Class`" = as.formula(~`UCA Class`*`Stress Class`*`Habitat Stability Class`),
    "`UCA Class`:`Habitat Stability Class`" = as.formula(~`UCA Class`*`Habitat Stability Class`),
    "`UCA Class`:`Stress Class`" = as.formula(~`UCA Class`*`Stress Class`),
    "`Habitat Stability Class`" = as.formula(~`Habitat Stability Class`),
    "`Stress Class`" = as.formula(~`Stress Class`),
    "`UCA Class`" = as.formula(~`UCA Class`),
    "`1`" = as.formula(~1)
  )
  
  fixed_list <- map(fixed_list,
                    ~{
                      emt95 <- emmeans::emtrends(mod,.x,var="Year_n",level=0.95) %>% # consider: addl.vars = "SiteCode" ,component="response"
                        data.frame(check.names = F) %>% 
                        mutate(level=0.95,term="Year_n") %>% 
                        rename(val=Year_n.trend)
                      
                      emt50 <- emmeans::emtrends(mod,.x,var="Year_n",level=0.5) %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.5,term="Year_n") %>% 
                        rename(val=Year_n.trend)
                      
                      emm95 <- emmeans::emmeans(mod,.x,level=0.95, at = list(Year_n = max(mod$frame$Year_n)),type="response") %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.95,term="(Intercept)") %>% 
                        rename(val=response)
                      emm50 <- emmeans::emmeans(mod,.x,level=0.50, at = list(Year_n = max(mod$frame$Year_n)),type="response") %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.5,term="(Intercept)") %>% 
                        rename(val=response)
                      
                      bind_rows(emm95,emm50,emt95,emt50) %>% 
                        rowwise() %>% 
                        mutate(grp=paste0(c_across(any_of(c("UCA Class","Stress Class","Habitat Stability Class"))),collapse = ":"))
                    })
  
  fixed_list <- bind_rows(fixed_list,.id="grpvar") %>% 
    mutate(component="fixed") 
  
  re_tbl <- as.data.frame(glmmTMB::ranef(mod, condVar = T)) %>% 
    filter(component == "cond") %>% 
    filter(grpvar == "`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode") %>% 
    mutate(SiteCode = str_split(grp,":",simplify = T)[,4]) %>% 
    rename(val=condval,
           SE=condsd)
  
  out_frame_raw <- distinct(mod$frame[,sapply(mod$frame,class)=="factor"]) %>% 
    mutate(val = NA_real_,
           SE = NA_real_)
  
  out_frame <- out_frame_raw
  
  for (i in unique(out_frame$SiteCode)) {
    sub_frame1 <- fixed_list %>% 
      filter(term == "Year_n") %>% 
      filter(level == 0.95) %>% 
      filter(
        `UCA Class` == out_frame[out_frame$SiteCode==i,"UCA Class"],
        `Habitat Stability Class` == out_frame[out_frame$SiteCode==i,"Habitat Stability Class"],
        `Stress Class` == out_frame[out_frame$SiteCode==i,"Stress Class"]
      )
    
    sub_frame2 <- re_tbl %>% 
      filter(SiteCode==i) %>% 
      filter(term == "Year_n")
    
    out_frame$val[out_frame$SiteCode==i] <- sum(sub_frame1$val,na.rm = T) + sum(sub_frame2$val,na.rm = T)
    out_frame$SE[out_frame$SiteCode==i] <- sqrt(sum(sub_frame1$SE^2,na.rm = T) + sum(sub_frame2$SE^2,na.rm = T))
    
  }
  
  out_frame <- out_frame %>% 
    mutate(
      pred_low = val - (qnorm(1 - 0.1 / 2) * SE),
      pred_high = val + (qnorm(1 - 0.1 / 2) * SE),
      pred_low_trend = val - (qnorm(1 - 0.5 / 2) * SE),
      pred_high_trend = val + (qnorm(1 - 0.5 / 2) * SE),
    ) %>% 
    mutate(
      term = "Year_n",
      component = "cond",
      grpvar = "`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode",
    )
  
  fixed_list_wide <- fixed_list %>% 
    pivot_wider(names_from = level, values_from = c(asymp.LCL,asymp.UCL)) %>% 
    rename(
      pred_low = asymp.LCL_0.95,
      pred_low_trend = asymp.LCL_0.5,
      pred_high = asymp.UCL_0.95,
      pred_high_trend = asymp.UCL_0.5
    )
  
  
  
  return(
    bind_rows(out_frame,fixed_list_wide) %>% 
      rename(value=val,
             group = grp,
             group_var = grpvar,
             sd = SE)
  )
  
}

summary_fun <- function(x){
  out <- map(
    c("Subshed","WSHED","UCA Class","Stress Class","UCA Class:Stress Class","UCA Class:Stress Class","UCA Class:Habitat Stability Class","UCA Class:Stress Class:Habitat Stability Class"),
    ~{x %>% 
        as_tibble() %>% 
        mutate(`UCA Class:Stress Class` = paste(`UCA Class`,`Stress Class`,sep = ".")) %>% 
        mutate(`UCA Class:Habitat Stability Class` = paste(`UCA Class`,`Habitat Stability Class`,sep = ".")) %>% 
        mutate(`UCA Class:Stress Class:Habitat Stability Class` = paste(`UCA Class`,`Stress Class`,`Habitat Stability Class`,sep = ".")) %>% 
        mutate(group=.x) %>% 
        rename_with(.cols=any_of(.x),~paste0("group_val")) %>% 
        mutate(group_val=as.character(group_val)) %>% 
        group_by(group,group_val) %>% 
        summarise(
          across(c(
            `Upstream Catchment Area`,
          ),
          ~mean(log10(.x),na.rm=T)^10
          ),
          across(c(
            `Reach Slope`,
          ),
          ~expm1(mean(log1p(.x),na.rm=T))
          ),
          across(c(
            `Baseflow Index`,
            `Upstream catchment developed (%)`,
            `Upstream catchment agricultural (%)`,
            `Upstream catchment stress index`
          ),
          ~median(.x,na.rm=T)
          ),
          .groups="drop"
        )
    }
  )
  
  out <- bind_rows(out)
  return(out)
}

# Function to extract model optim

glmmTMB_get_optimum <- function(mod) {
  # Fixed effects
  starting_point <- list(
    beta = glmmTMB::fixef(mod)$cond,
    betazi = glmmTMB::fixef(mod)$zi,
    betadisp = glmmTMB::fixef(mod)$disp
  )
  #Random effects
  ref <- glmmTMB::ranef(mod)
  if (length(ref$cond) > 0) {
    starting_point$theta = ref$cond
  }
  if (length(ref$zi) > 0) {
    starting_point$thetazi = ref$zi
  }
  if (length(ref$disp) > 0) {
    starting_point$thetadisp = ref$disp
  }
  return(starting_point)
}

# Function to generate heatmaps
heatmap_fn<-function(x_df,
                     y_df=NULL,
                     alpha=0.1,
                     sig_est=0.6,
                     method="spearman",
                     scale_grad=ggplot2::scale_fill_gradient2(
                       low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                       mid = "white",
                       high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                       limits=c(-1,1)
                     ),
                     ...){
  
  if (is.null(y_df)){
    y_df<-x_df
  }
  cor_list<-lapply(x_df,function(x) lapply(y_df,function(y) broom::tidy(cor.test(x,y,method=method,...))))
  
  pcorr_out<-purrr::map_dfr(cor_list,.id="x_var",~purrr::map_dfr(.,dplyr::bind_rows,.id="y_var"))
  
  pcorr_out <- pcorr_out %>%
    mutate(p.value = case_when(
      estimate >= 1 - .Machine$double.eps ~ NA,
      estimate <= -1 + .Machine$double.eps~ NA,
      T ~ p.value
    )) %>%
    mutate(estimate = case_when(
      estimate >= 1 - .Machine$double.eps ~ NA,
      estimate <= -1 + .Machine$double.eps~ NA,
      T ~ estimate
    )) 
  
  rho<-pcorr_out %>%
    select(x_var,y_var,estimate) %>%
    spread(x_var,estimate) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var") %>%
    textshape::cluster_matrix(method ="ward.D2")
  
  pv<-pcorr_out %>%
    mutate(p.value=case_when(
      is.na(p.value) ~ "darkgray",
      p.value>alpha ~ "darkgray",
      p.value<=alpha & abs(estimate)>=sig_est ~ "darkred",
      p.value<=alpha ~ "orange",
      T ~ "darkgray"
    )) %>%
    select(x_var,y_var,p.value) %>%
    spread(x_var,p.value) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var")
  
  pv<-pv[rownames(rho),colnames(rho)]
  
  plot_out <- pcorr_out %>%
    select(x_var,y_var,estimate) %>% 
    spread(x_var,estimate) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var") %>%
    heatmaply::ggheatmap(hclust_method ="ward.D2",
                         seriate="GW",
                         margins =c(40,40,40,40),
                         branches_lwd = 0.25,
                         grid_gap=3,
                         heatmap_layers=list(
                           theme(
                             axis.text.y = element_text(size = 8)   ,
                             axis.text.x = element_text(
                               size = 8,
                               hjust = 1,
                               vjust = 1,
                               angle = 45
                             ),
                             legend.position = "none"),
                           geom_tile(data=pv %>%
                                       rownames_to_column("Endpoint") %>%
                                       gather("Parameter","Value",-Endpoint),
                                     aes(x=Parameter,y=Endpoint),
                                     #alpha=0,
                                     colour=pv %>%
                                       rownames_to_column("Endpoint") %>%
                                       gather("Parameter","Value",-Endpoint) %>%
                                       pull(Value),
                                     width=0.975, height=0.975,
                                     size=0.3,
                                     fill=NA,
                                     inherit.aes = F),
                           guides(colour="none")
                         ),
                         scale_fill_gradient_fun = scale_grad)
  
  cowplot::ggdraw(plot_out)
  
}

# transformer for car::logit
clogit <- scales::new_transform(
  name="clogit",
  transform = function(x) {
    x[x<0] <- 0
    x[x>1] <- 1
    car::logit(x,T,0.025)
  },
  inverse = function(x) {
    a <- 0.025
    a <- (1-2*a)
    a <- (a*(1+exp(x))+(exp(x)-1))/(2*a*(1+exp(x)))
    zapsmall(a)
  },
  domain = c(0, 1)
)

# Functions to plot axis on pseudo-log scales
force_all <- function(...) list(...) # Helper function defined in scales

pseudo_log <- function (val, sigma = 1, base = exp(1)) {
  asinh(val/(2*sigma))/log(base)
}

pseudo_log_breaks <- function(n = 5, sigma = 1, base = 10) {
  force_all(n, base)
  n_default = n
  
  force_all <- function(...) list(...) # Helper function defined in scales
  
  pseudo_log <- function (val, sigma = 1, base = exp(1)) {
    asinh(val/(2*sigma))/log(base)
  }
  
  function(x, n = n_default) {
    rng <- pseudo_log(range(x, na.rm = TRUE), sigma = sigma, base = base)
    min <- floor(rng[1])
    max <- ceiling(rng[2])
    
    if (max == min) return(base^min)
    
    by <- floor((max - min) / n) + 1
    breaks <- base^seq(min, max, by = by)
    relevant_breaks <- base^rng[1] <= breaks & breaks <= base^rng[2]
    if (sum(relevant_breaks) >= (n - 2)) return(breaks)
    
    # the easy solution to get more breaks is to decrease 'by'
    while (by > 1) {
      by <- by - 1
      breaks <- base^seq(min, max, by = by)
      relevant_breaks <- base^rng[1] <= breaks & breaks <= base^rng[2]
      if (sum(relevant_breaks) >= (n - 2)) return(breaks)
    }
    scales:::log_sub_breaks(rng, n = n, base = base)
  }
}

# Helper function to extract a legend
get_legend <- function(p, position = NULL){
  
  if(.is_list(p)){
    continue <- TRUE
    i <- 1
    while(i <= length(p) & continue){
      leg <- .get_legend(p[[i]], position = position)
      if(!is.null(leg)) continue <- FALSE
      i <- i+1
    }
  }
  else{
    leg <- .get_legend(p, position = position)
  }
  leg
}

.is_list <- function(x){
  inherits(x, "list")
}


# Return legend for one plot
.get_legend <- function(p, position = NULL){
  
  if(is.null(p)) return(NULL)
  if(!is.null(position)){
    p <- p + theme(legend.position = position)
  }
  tmp <- ggplot_gtable(ggplot_build(p))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  if(length(leg) > 0) leg <- tmp$grobs[[leg]]
  else leg <- NULL
  leg
}

col2hex <- function(color, alpha = 1) {
  
  if (alpha < 0 | alpha > 1)
    stop('alpha must be between 0 and 1.')
  
  if (missing(color))
    stop('Need a color')
  
  if (!is.character(color))
    stop(strwrap('Elements must be character string as a named R color or
         hex (e.g. "#ffffff")'))
  
  if (!color %in% colors()) {
    if (!nchar(color) == 7 | !grepl('^#', color))
      stop(strwrap('color must be an R color, e.g. see colors(),
                   or a hex of the form #ff5500'))
  }
  
  rgb(t(col2rgb(color)/255), alpha = alpha)
}

# Calculate fish endpoints
calc_fish_ep<-function(df) {
  
  ep_tbl<-df %>% 
    group_by(SampleEventID,SampleDate) %>% 
    mutate(
      # General Community
      Comm_Biomass=sum(EstimatedBiomass,na.rm=T),
      Comm_Abundance=sum(NumberOfFish,na.rm=T),
      Comm_Biomass_SATI=sum(EstimatedBiomass[!is.na(OPT_TEMP)],na.rm=T),
      Comm_Abundance_SATI=sum(NumberOfFish[!is.na(OPT_TEMP)],na.rm=T),
      Comm_Richness=length(unique(SpeciesCode)),
      EstimatedBiomass_perc=EstimatedBiomass/Comm_Biomass,
      NumberOfFish_perc=NumberOfFish/Comm_Abundance,
      EstimatedBiomass_perc_SATI=EstimatedBiomass/Comm_Biomass_SATI,
      NumberOfFish_perc_SATI=NumberOfFish/Comm_Abundance_SATI,
      SATI_EstimatedBiomass_perc=EstimatedBiomass_perc_SATI*OPT_TEMP,
      SATI_NumberOfFish_perc=NumberOfFish_perc_SATI*OPT_TEMP,
    ) %>% 
    summarise(
      # General Community
      Comm_Biomass=head(Comm_Biomass,1),
      Comm_Abundance=head(Comm_Abundance,1),
      Comm_Richness=head(Comm_Richness,1),
      # Thermal Regime 
      Therm_SATI_Bioperc=sum(SATI_EstimatedBiomass_perc,na.rm=T),
      Therm_SATI_Abuperc=sum(SATI_NumberOfFish_perc,na.rm=T),
      
      Therm_cold_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime=="coldwater"],na.rm=T),
      Therm_cool_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime=="coolwater"],na.rm=T),
      Therm_coolcold_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T),
      Therm_warm_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("warmwater")],na.rm=T),
      Therm_cold_Abuperc=sum(NumberOfFish_perc[Thermal.Regime=="coldwater"],na.rm=T),
      Therm_cool_Abuperc=sum(NumberOfFish_perc[Thermal.Regime=="coolwater"],na.rm=T),
      Therm_coolcold_Abuperc=sum(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T),
      Therm_warm_Abuperc=sum(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")],na.rm=T),
      Therm_cold_Bio=sum(EstimatedBiomass_perc[Thermal.Regime=="coldwater"],na.rm=T)*Comm_Biomass,
      Therm_cool_Bio=sum(EstimatedBiomass_perc[Thermal.Regime=="coolwater"],na.rm=T)*Comm_Biomass,
      Therm_coolcold_Bio=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T)*Comm_Biomass,
      Therm_warm_Bio=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("warmwater")],na.rm=T)*Comm_Biomass,
      Therm_cold_Abu=sum(NumberOfFish_perc[Thermal.Regime=="coldwater"],na.rm=T)*Comm_Abundance,
      Therm_cool_Abu=sum(NumberOfFish_perc[Thermal.Regime=="coolwater"],na.rm=T)*Comm_Abundance,
      Therm_coolcold_Abu=sum(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T)*Comm_Abundance,
      Therm_warm_Abu=sum(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")],na.rm=T)*Comm_Abundance,
      Therm_warm_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")]),
      Therm_cold_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coldwater")]),
      Therm_cool_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coolwater")]),
      Therm_coolcold_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")]),
      # Trophic Guild 
      Troph_invert_Bioperc=sum(EstimatedBiomass_perc[grepl("invertivore",Trophic.Class)],na.rm=T),
      Troph_detrit_Bioperc=sum(EstimatedBiomass_perc[grepl("detritivore",Trophic.Class)],na.rm=T),
      Troph_carn_Bioperc=sum(EstimatedBiomass_perc[grepl("carnivore",Trophic.Class)],na.rm=T),
      Troph_herb_Bioperc=sum(EstimatedBiomass_perc[grepl("herbivore",Trophic.Class)],na.rm=T),
      Troph_plank_Bioperc=sum(EstimatedBiomass_perc[grepl("planktivore",Trophic.Class)],na.rm=T),
      Troph_invert_Bio=sum(EstimatedBiomass_perc[grepl("invertivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_detrit_Bio=sum(EstimatedBiomass_perc[grepl("detritivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_carn_Bio=sum(EstimatedBiomass_perc[grepl("carnivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_herb_Bio=sum(EstimatedBiomass_perc[grepl("herbivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_plank_Bio=sum(EstimatedBiomass_perc[grepl("planktivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_invert_Abuperc=sum(NumberOfFish_perc[grepl("invertivore",Trophic.Class)],na.rm=T),
      Troph_detrit_Abuperc=sum(NumberOfFish_perc[grepl("detritivore",Trophic.Class)],na.rm=T),
      Troph_carn_Abuperc=sum(NumberOfFish_perc[grepl("carnivore",Trophic.Class)],na.rm=T),
      Troph_herb_Abuperc=sum(NumberOfFish_perc[grepl("herbivore",Trophic.Class)],na.rm=T),
      Troph_plank_Abuperc=sum(NumberOfFish_perc[grepl("planktivore",Trophic.Class)],na.rm=T),
      Troph_invert_Abu=sum(NumberOfFish_perc[grepl("invertivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_detrit_Abu=sum(NumberOfFish_perc[grepl("detritivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_carn_Abu=sum(NumberOfFish_perc[grepl("carnivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_herb_Abu=sum(NumberOfFish_perc[grepl("herbivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_plank_Abu=sum(NumberOfFish_perc[grepl("planktivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_invert_Rich=length(NumberOfFish_perc[grepl("invertivore",Trophic.Class)]),
      Troph_detrit_Rich=length(NumberOfFish_perc[grepl("detritivore",Trophic.Class)]),
      Troph_carn_Rich=length(NumberOfFish_perc[grepl("carnivore",Trophic.Class)]),
      Troph_herb_Rich=length(NumberOfFish_perc[grepl("herbivore",Trophic.Class)]),
      Troph_plank_Rich=length(NumberOfFish_perc[grepl("planktivore",Trophic.Class)]),
      # Tolerance
      Tol_tol_Bioperc=sum(EstimatedBiomass_perc[Tolerance=="tolerant"],na.rm=T),
      Tol_intol_Bioperc=sum(EstimatedBiomass_perc[Tolerance=="intolerant"],na.rm=T),
      Tol_tol_Abuperc=sum(NumberOfFish_perc[Tolerance=="tolerant"],na.rm=T),
      Tol_intol_Abuperc=sum(NumberOfFish_perc[Tolerance=="intolerant"],na.rm=T),
      Tol_tol_Bio=sum(EstimatedBiomass_perc[Tolerance=="tolerant"],na.rm=T)*Comm_Biomass,
      Tol_intol_Bio=sum(EstimatedBiomass_perc[Tolerance=="intolerant"],na.rm=T)*Comm_Biomass,
      Tol_tol_Abu=sum(NumberOfFish_perc[Tolerance=="tolerant"],na.rm=T)*Comm_Abundance,
      Tol_intol_Abu=sum(NumberOfFish_perc[Tolerance=="intolerant"],na.rm=T)*Comm_Abundance,
      Tol_tol_Rich=length(NumberOfFish_perc[Tolerance=="tolerant"]),
      Tol_intol_Rich=length(NumberOfFish_perc[Tolerance=="intolerant"]),
      # Individual Species - 
      Spec_RBD_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==337],na.rm=T),
      Spec_RBD_Abuperc=sum(NumberOfFish_perc[SpeciesCode==337],na.rm=T),
      Spec_BKT_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==80],na.rm=T),
      Spec_BKT_Abuperc=sum(NumberOfFish_perc[SpeciesCode==80],na.rm=T),
      Spec_GDF_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==181],na.rm=T),
      Spec_GDF_Abuperc=sum(NumberOfFish_perc[SpeciesCode==181],na.rm=T),
      Spec_CMC_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==186],na.rm=T),
      Spec_CMC_Abuperc=sum(NumberOfFish_perc[SpeciesCode==186],na.rm=T),
      Spec_RDG_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==366],na.rm=T),
      Spec_RDG_Abuperc=sum(NumberOfFish_perc[SpeciesCode==366],na.rm=T),
      Spec_CKC_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==212],na.rm=T),
      Spec_CKC_Abuperc=sum(NumberOfFish_perc[SpeciesCode==212],na.rm=T),
      Spec_SKP_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(381,382)],na.rm=T),
      Spec_SKP_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(381,382)],na.rm=T),
      Spec_LMP_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(11:14)],na.rm=T),
      Spec_LMP_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(11:14)],na.rm=T),
      Spec_DRT_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(335:340,342)],na.rm=T),
      Spec_DRT_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(335:340,342)],na.rm=T),
      Spec_CEB_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T),
      Spec_CEB_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T),
      Spec_CTF_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(231:233)],na.rm=T),
      Spec_CTF_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(231:233)],na.rm=T),
      Spec_SAL_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(71:83)],na.rm=T),
      Spec_SAL_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(71:83)],na.rm=T),
      
      Spec_RBD_Bio=sum(EstimatedBiomass_perc[SpeciesCode==337],na.rm=T)*Comm_Biomass,
      Spec_RBD_Abu=sum(NumberOfFish_perc[SpeciesCode==337],na.rm=T)*Comm_Abundance,
      Spec_BKT_Bio=sum(EstimatedBiomass_perc[SpeciesCode==80],na.rm=T)*Comm_Biomass,
      Spec_BKT_Abu=sum(NumberOfFish_perc[SpeciesCode==80],na.rm=T)*Comm_Abundance,
      Spec_GDF_Bio=sum(EstimatedBiomass_perc[SpeciesCode==181],na.rm=T)*Comm_Biomass,
      Spec_GDF_Abu=sum(NumberOfFish_perc[SpeciesCode==181],na.rm=T)*Comm_Abundance,
      Spec_CMC_Bio=sum(EstimatedBiomass_perc[SpeciesCode==186],na.rm=T)*Comm_Biomass,
      Spec_CMC_Abu=sum(NumberOfFish_perc[SpeciesCode==186],na.rm=T)*Comm_Abundance,
      Spec_RDG_Bio=sum(EstimatedBiomass_perc[SpeciesCode==366],na.rm=T)*Comm_Biomass,
      Spec_RDG_Abu=sum(NumberOfFish_perc[SpeciesCode==366],na.rm=T)*Comm_Abundance,
      Spec_CKC_Bio=sum(EstimatedBiomass_perc[SpeciesCode==212],na.rm=T)*Comm_Biomass,
      Spec_CKC_Abu=sum(NumberOfFish_perc[SpeciesCode==212],na.rm=T)*Comm_Abundance,
      Spec_SKP_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(381,382)],na.rm=T)*Comm_Biomass,
      Spec_SKP_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(381,382)],na.rm=T)*Comm_Abundance,
      Spec_LMP_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(11:14)],na.rm=T)*Comm_Biomass,
      Spec_LMP_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(11:14)],na.rm=T)*Comm_Abundance,
      Spec_DRT_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(335:340,342)],na.rm=T)*Comm_Biomass,
      Spec_DRT_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(335:340,342)],na.rm=T)*Comm_Abundance,
      Spec_CEB_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T)*Comm_Biomass,
      Spec_CEB_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T)*Comm_Abundance,
      Spec_CTF_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(231:233)],na.rm=T)*Comm_Biomass,
      Spec_CTF_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(231:233)],na.rm=T)*Comm_Abundance,
      Spec_SAL_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(71:83)],na.rm=T)*Comm_Biomass,
      Spec_SAL_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(71:83)],na.rm=T)*Comm_Abundance,
      # Environment
      Env_benthic_Bioperc=sum(EstimatedBiomass_perc[Environment=="benthic"],na.rm=T),
      Env_pelagic_Bioperc=sum(EstimatedBiomass_perc[Environment=="pelagic"],na.rm=T),
      Env_benthopelagic_Bioperc=sum(EstimatedBiomass_perc[Environment=="benthopelagic"],na.rm=T),
      Env_benthic_Abuperc=sum(NumberOfFish_perc[Environment=="benthic"],na.rm=T),
      Env_pelagic_Abuperc=sum(NumberOfFish_perc[Environment=="pelagic"],na.rm=T),
      Env_benthopelagic_Abuperc=sum(NumberOfFish_perc[Environment=="benthopelagic"],na.rm=T),
      Env_benthic_Bio=sum(EstimatedBiomass_perc[Environment=="benthic"],na.rm=T)*Comm_Biomass,
      Env_pelagic_Bio=sum(EstimatedBiomass_perc[Environment=="pelagic"],na.rm=T)*Comm_Biomass,
      Env_benthopelagic_Bio=sum(EstimatedBiomass_perc[Environment=="benthopelagic"],na.rm=T)*Comm_Biomass,
      Env_benthic_Abu=sum(NumberOfFish_perc[Environment=="benthic"],na.rm=T)*Comm_Abundance,
      Env_pelagic_Abu=sum(NumberOfFish_perc[Environment=="pelagic"],na.rm=T)*Comm_Abundance,
      Env_benthopelagic_Abu=sum(NumberOfFish_perc[Environment=="benthopelagic"],na.rm=T)*Comm_Abundance,
      Env_benthic_Rich=length(NumberOfFish_perc[Environment=="benthic"]),
      Env_pelagic_Rich=length(NumberOfFish_perc[Environment=="pelagic"]),
      Env_benthopelagic_Rich=length(NumberOfFish_perc[Environment=="benthopelagic"]),
      # Reproductive Guild
      ## Guarders
      Rep_guard_Bioperc=sum(EstimatedBiomass_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T),
      Rep_guard_Abuperc=sum(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T),
      Rep_guard_Bio=sum(EstimatedBiomass_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_guard_Abu=sum(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_guard_Rich=length(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)]),
      ## Nesting?
      Rep_OSS_Bioperc=sum(EstimatedBiomass_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T),
      Rep_OSS_Abuperc=sum(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T),
      Rep_OSS_Bio=sum(EstimatedBiomass_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_OSS_Abu=sum(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_OSS_Rich=length(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)]),
      Rep_BrHid_Bioperc=sum(EstimatedBiomass_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T),
      Rep_BrHid_Abuperc=sum(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T),
      Rep_BrHid_Bio=sum(EstimatedBiomass_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_BrHid_Abu=sum(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_BrHid_Rich=length(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)]),
      ## Lithophilic
      Rep_Lith_Bioperc=sum(EstimatedBiomass_perc[grepl("Litho",Reproductive.Guild)],na.rm=T),
      Rep_Lith_Abuperc=sum(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)],na.rm=T),
      Rep_Lith_Bio=sum(EstimatedBiomass_perc[grepl("Litho",Reproductive.Guild)],na.rm=T)*Comm_Biomass,
      Rep_Lith_Abu=sum(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)],na.rm=T)*Comm_Abundance,
      Rep_Lith_Rich=length(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)]),
      .groups="drop"
    ) 
  
  return(ep_tbl)
}

2 Download Data

source("00_Functions/FWIS_Functions.R")
library(tidyverse)

username<-Sys.getenv("FWIS_username")
password<-Sys.getenv("FWIS_password")

fwis_login<-wideLogin("http://www.comap.ca/fwis/wideR.php",username,password)


# Habitat Data ------------------------------------------------------------

hab_tbls<-c("tblChannelStructureSummary")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",hab_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=hab_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/tblChannelStructureSummary.csv")

# Fish Data ---------------------------------------------------------------

fish_tbls<-c("tblFishSummaryOfTotalCatches")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",fish_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=fish_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/tblFishSummaryOfTotalCatches.csv")

# Fish Habitat Data -------------------------------------------------------

fish_tbls<-c("lkpSpeciesHabitatRequirement")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",fish_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=fish_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/lkpSpeciesHabitatRequirement.csv")

# Benthic Invertebrate Data -----------------------------------------------

# Benthic Invertebrate Download -------------------------------------------

tables <- unique(wideDataSelect(1,"tblDataAPI", list("tbldUserTblName"), orderBy = list("tbldUserTblName")))

bic_tbls<-tables$tbldUserTblName[grepl("benthic",tolower(tables$tbldUserTblName)) & grepl("tbl",tolower(tables$tbldUserTblName))]

bic_all<-lapply(setNames(bic_tbls,bic_tbls),function(x){
  columns <- wideDataSelect(1,
                            "tblDataAPI",
                            list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                            where=list(tbldUserTblName=c("=",x)),
                            orderBy=list("tbldOrdering"))
  
  wideDataSelect( 1, table=x, colsSelected=as.list(columns$tbldColName) )
  
})

saveRDS(bic_all,"data/BIC_ALL.rds")

obbn_tbls<-tables$tbldUserTblName[grepl("obbn",tolower(tables$tbldUserTblName)) & grepl("tbl",tolower(tables$tbldUserTblName))]

obbn_all<-lapply(setNames(obbn_tbls,obbn_tbls),function(x){
  columns <- wideDataSelect(1,
                            "tblDataAPI",
                            list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                            where=list(tbldUserTblName=c("=",x)),
                            orderBy=list("tbldOrdering"))
  
  columns$tbldColName <- gsub("Â", "", columns$tbldColName)
  
  wideDataSelect( 1, table=x, colsSelected=as.list(columns$tbldColName) )
  
})

saveRDS(obbn_all,"data/OBBN_ALL.rds")

3 Format Raw Data

source("00_source.R")
library(openxlsx)


# Read in Data ------------------------------------------------------------
dat <- read_csv("data/tblChannelStructureSummary.csv", locale = locale(encoding = "latin1"))


# Add additional endpoints ------------------------------------------------
dat <- dat %>% 
  mutate(
    TotalCoverTypePCWood = EmbCoverTypePCWood + UnembCoverTypePCWood,
    TotalCoverTypePCRndRock = EmbCoverTypePCRndRock + UnembCoverTypePCRndRock,
    TotalCoverTypePCFltRock = EmbCoverTypePCFltRock + UnembCoverTypePCFltRock,
    TotalCoverTypePCMacrophyte = EmbCoverTypePCMacrophyte + UnembCoverTypePCMacrophyte,
    TotalCoverTypePCBank = EmbCoverTypePCBank + UnembCoverTypePCBank,
  ) %>% 
  mutate(
    TotalVegTypePCFL = EmbVegTypePCFL + UnembVegTypePCFL,
    TotalVegTypePCAL = EmbVegTypePCAL + UnembVegTypePCAL,
    TotalVegTypePCSS = EmbVegTypePCSS + UnembVegTypePCSS,
    TotalVegTypePCMC = EmbVegTypePCMC + UnembVegTypePCMC,
  )

# Exclude test data
dat <- dat %>% 
  select(all_of(c(unlist(ID_cols),unlist(met_cols),StreamType = "substrateQualityRatingCode"))) %>% 
  filter(!StreamCode %in% c("TEST")) %>% 
  filter(!grepl("Training",OrgNameList)) %>% 
  filter(SiteCode != "NCD21-1C") #%>% # TRCA site not within bounds of TRCA shape file
  #filter(SiteCode != "PT003WM") %>%  # incorrect coordinates
  #filter(SiteCode != "PT002WM") 
# filter(month(SampleDate) != "1",
#        month(SampleDate) != "2",
#        month(SampleDate) != "12")

# Remove samples without location information
dat <- dat  %>% 
  rename_with(~gsub("^.*?\\.","",.x)) %>% 
  mutate(Longitude = as.numeric(Longitude),
         Latitude = as.numeric(Latitude)) %>% 
  group_by(SiteCode) %>% 
  mutate(Longitude = median(Longitude,na.rm = T),
         Latitude = median(Latitude,na.rm = T),
  ) %>% 
  filter(!is.na(Longitude),
         !is.na(Latitude)) %>% 
  ungroup()

# Add utility columns
dat <- dat %>% 
  group_by(!!!syms(ID_cols[!ID_cols %in% c("ID","SampleDate","Sample","Longitude","Latitude")]))  %>% 
  mutate(
    Month = month(SampleDate),
    Year = year(SampleDate),
    `Sample Season` = case_when(
      Month %in% 3:5 ~ "Spring",
      Month %in% 6:8 ~ "Summer",
      Month %in% 9:11 ~ "Fall",
      T ~ "Winter"
    ),
    `Sample Years` = paste(sort(unique(Year)),collapse = ", "),
    `Sample Months` = paste(sort(unique(Month)),collapse = ", "),
    `Years of Data` = length(unique(Year)),
    `Year-Month`=format(SampleDate,"%Y-%m"),
    `Sample Year-Months` = paste(sort(unique(`Year-Month`)),collapse = ", "),
    `Sample Year-Season` = paste(sort(unique(`Year-Month`)),collapse = ", "),
    `TRCA` = OrgNameList == "Toronto and Region Conservation Authority",
    `Sample Season`=factor(`Sample Season`,levels=c("Sping","Summer","Fall","Winter"))
  ) %>% 
  filter(Year < 2025) %>% 
  ungroup() %>% 
  group_by(!!!syms(c(ID_cols[!ID_cols %in% c("ID","SampleDate","Sample","Longitude","Latitude")],"Year"))) %>% 
  mutate( 
    `Sampled Days within Year` = length(unique(SampleDate)),
    `Sampled Seasons within Year` = length(unique(`Sample Season`)),
  ) %>% 
  ungroup() %>% 
  rowwise() %>% 
  mutate(
    `Number of Missing Endpoints` = sum(is.na(c_across(unlist(sapply(met_cols,names),use.names = F))))
  ) %>% 
  ungroup() %>% 
  mutate(
    Year_n = (max(Year)-min(Year))-(max(Year)-Year)+1,
    year_class = case_when(
      Year %in% 1995:2000 ~ "1995-2000",
      Year %in% 2001:2006 ~ "2001-2006",
      Year %in% 2007:2012 ~ "2007-2012",
      Year %in% 2013:2018 ~ "2013-2018",
      Year %in% 2019:2024 ~ "2019-2024",
    )
  ) %>% 
  mutate(
    year_class = factor(year_class,levels=rev(c("1995-2000","2001-2006","2007-2012","2013-2018","2019-2024")))
  )


# Duplicate TRCA data in full analysis subset -----------------------------
dat <- bind_rows(
  dat %>% mutate(TRCA=F),
  dat %>% filter(TRCA) %>% mutate(TRCA=T)
)


# Filter Data -------------------------------------------------------------
dat <- dat %>% 
  filter(`Years of Data` > 2 & length(unique(year_class)) > 1) %>% 
  group_by(SiteCode,TRCA) %>% 
  filter(
    !TRCA | `Years of Data` > 3,
    !TRCA | any(Year > 2018),
    !TRCA | length(unique(year_class)) > 1
  ) %>% 
  ungroup() 



if (F){
  # Get Ranges -------------------------------------------------------------
  
  dat %>% 
    pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
                 names_to = "Endpoint") %>% 
    ggplot(aes(x=value)) +
    geom_histogram() +
    facet_wrap(~Endpoint,scales="free")
  
  dat %>% 
    pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
                 names_to = "Endpoint") %>% 
    group_by(Endpoint) %>% 
    summarise(
      min = min(value,na.rm = T),
      max = max(value,na.rm = T)
    ) %>% 
    print(n=Inf)
  
}

# Convert Percent to Proportion -------------------------------------------

dat <- dat %>% 
  mutate(
    across(
      contains(c("Cover","Percent")),
      ~ .x/100
    )
  )

# Replace problematic values ------------------------------------------------


dat <- dat %>% 
  mutate(
    `Habitat Stability Score` = `Habitat Stability Score` - 0.5,
    `Habitat Homogeneity` = 1-`Habitat Homogeneity`,
    across(
      c(D15,D50,D85,D50Max),
      ~ case_when(
        .x < 0 ~ NA_real_,
        .x <= 2 ~ 0,
        .x >= 1000 ~ NA_real_,
        T ~ .x
      )
    ),
    D15 = case_when(
      D15 > 50 ~ NA_real_,
      T ~ D15
    ),
    across(
      c(`Sediment Transport Potential`),
      ~ case_when(
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(`Sediment Sorting Index`),
      ~ case_when(
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(
        contains(c("Cover","Percent")),
        `Habitat Stability Score`,
        `Habitat Homogeneity`,
        `Bank Stability - Mean`
      ),
      ~ case_when(
        .x >= 1 ~ quantile(.x[.x<1 & .x>0],0.99,na.rm=T),
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(
        contains(c("Cover","Percent"))
      ),
      ~ case_when(
        .x < 0.02 ~ 0,
        T ~ .x
      )
    )
  )


# Convert to sf, and attribute with watershed info

dat <- dat %>% 
  sf::st_as_sf(coords=c("Longitude","Latitude"),crs="WGS84",remove = F)


# Geospatial analysis -----------------------------------------------------

TRCA_subwats <- sf::read_sf("data/Subwatersheds_TRCA.geojson")

sf::sf_use_s2(F) # there is an invalid geometry

unzip("data/AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package06_LakeSuperior.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers.zip",exdir=tempdir())
strm_lns <- list(
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w03_Lake_Ontario_West/w03_AEC_V3_Core_Rev0.gpkg"),
    "w03_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w06_Lake_Ontario_East/w06_AEC_V3_Core_Rev0.gpkg"),
    "w06_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w04_Lake_Ontario_Central/w04_AEC_V3_Core_Rev0.gpkg"),
    "w04_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package06_LakeSuperior/w23_Lake_Superior_North_West/w23_AEC_V3_Core_Rev0.gpkg"),
    "w23_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers/w08_Ottawa_River_Lower/w08_AEC_V3_Core_Rev0.gpkg"),
    "w08_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers/w09_Ottawa_River_Central_South/w09_AEC_V3_Core_Rev0.gpkg"),
    "w09_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth/w14_Georgian_Bay_South_Simcoe/w14_AEC_V3_Core_Rev0.gpkg"),
    "w14_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth/w15_Georgian_Bay_Central/w15_AEC_V3_Core_Rev0.gpkg"),
    "w15_Reach_SimpleGeometry_Core")
)

strm_lns <- bind_rows(strm_lns)

strm_lns <- strm_lns %>% 
  filter(Network_Line_Type != "Shoreline Virtual Connector",
         Network_Line_Type != "Virtual Connector") %>% 
  sf::st_transform(sf::st_crs(dat))

dat <- sf::st_join(dat,
                   TRCA_subwats %>% select(WSHED,Subshed)
)

dat <- sf::st_join(dat,
                   strm_lns %>% select(WorkUnitID,ProvReachID,`Strahler Order`=Strahler_Order,`Upstream Catchment Area`=Upstream_Catchment_Area,`Baseflow Index`=BFI_RCA,`Reach Slope`=Slope_ReachChannel_Percent,Turbidity_percUpstreamChannel_T,Turbidity,Flow_Velocity),
                   join =nngeo::st_nn,k=1,maxdist=Inf,parallel=8
)

dat <- dat %>% 
  left_join(
    read_csv("data/AEC_supplementary_MWP_human_disturbance_gradients.csv") %>% 
      select(ProvReachID = ReachID, `Upstream catchment developed (%)`,
             `Upstream catchment agricultural (%)`,
             `Upstream catchment stress index`,
             `Upstream catchment disturbance status`)
  )

dat <- dat %>% 
  mutate(
    WSHED = case_when(
      is.na(WSHED) ~ as.character(WorkUnitID),
      T ~ WSHED
    ),
    Subshed = case_when(
      is.na(Subshed) ~ OrgNameList,
      T ~ Subshed
    )
  )

if (F) mapview::mapview(list(strm_lns,dat))


# Summarize TRCA Streamlines ----------------------------------------------


TRCA_strm_lns <- strm_lns %>% 
  filter(WorkUnitID %in% unique(dat$WorkUnitID)) %>% 
  sf::st_join(TRCA_subwats %>% select(WSHED,Subshed)) %>% 
  left_join(
    read_csv("data/AEC_supplementary_MWP_human_disturbance_gradients.csv") %>% 
      select(ProvReachID = ReachID, `Upstream catchment developed (%)`,
             `Upstream catchment agricultural (%)`,
             `Upstream catchment stress index`,
             `Upstream catchment disturbance status`)
  ) %>% 
  rename(
    `Strahler Order` = Strahler_Order,
    `Upstream Catchment Area`=Upstream_Catchment_Area,
    `Baseflow Index`=BFI_RCA,
    `Reach Slope`=Slope_ReachChannel_Percent
  )

TRCA_strm_lns <- bind_rows(
  TRCA_strm_lns %>% 
    filter(!is.na(Subshed)) %>% 
    mutate(TRCA=T),
  TRCA_strm_lns %>% 
    mutate(TRCA=F) %>% 
    mutate(WSHED = as.character(WorkUnitID),
           Subshed = as.character(WorkUnitID))
)

TRCA_strm_lns <- TRCA_strm_lns %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

TRCA_strm_lns <- TRCA_strm_lns %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

levels(TRCA_strm_lns$`Stress Class`) <- c("Low","Medium","High")
TRCA_strm_lns$`Stress Class` <- as.character(TRCA_strm_lns$`Stress Class`)
TRCA_strm_lns$`Stress Class`[is.na(TRCA_strm_lns$`Stress Class`)] <- "North"
TRCA_strm_lns$`Stress Class` <- factor(TRCA_strm_lns$`Stress Class`,levels = c("Low","Medium","High","North"))

TRCA_strm_lns <- map(
  c("Subshed","WSHED","Strahler Order","Stress Class","Strahler Order:Stress Class"),
  ~{TRCA_strm_lns %>% 
      as_tibble() %>% 
      mutate(`Strahler Order:Stress Class` = paste(`Strahler Order`,`Stress Class`,sep = ".")) %>% 
      mutate(group=.x) %>% 
      rename_with(.cols=any_of(.x),~paste0("group_val")) %>% 
      mutate(group_val=as.character(group_val)) %>% 
      group_by(TRCA,group,group_val) %>% 
      summarise(
        across(c(
          `Upstream Catchment Area`,
        ),
        ~mean(log10(.x),na.rm=T)^10
        ),
        across(c(
          `Reach Slope`,
        ),
        ~expm1(mean(log1p(.x),na.rm=T))
        ),
        across(c(
          `Baseflow Index`,
          `Upstream catchment developed (%)`,
          `Upstream catchment agricultural (%)`,
          `Upstream catchment stress index`
        ),
        ~median(.x,na.rm=T)
        )
      )
  }
)

TRCA_strm_lns <- bind_rows(TRCA_strm_lns)

write_csv(TRCA_strm_lns,"data/predictor_summaries.csv")

# Convert locations to factors --------------------------------------------

dat <- dat %>% 
  group_by(WSHED) %>% 
  mutate(Longitude_w = median(Longitude),
         Latitude_w = median(Latitude)) %>% 
  group_by(Subshed) %>% 
  mutate(Longitude_s = median(Longitude),
         Latitude_s = median(Latitude)) %>% 
  mutate(
    WSHED = factor(WSHED,levels=unique(WSHED[order(Longitude_w)])),
    Subshed = factor(Subshed,levels=unique(Subshed[order(Longitude_s)])),
    SiteCode = factor(SiteCode,levels=unique(SiteCode[order(Longitude_s)]))
  ) %>% 
  ungroup()


# Bin  streams by cross-sectional size ------------------------------------

dat <- dat %>% 
  mutate(
    Area = `Width - Mean`*`Depth - Mean`
  ) %>% 
  group_by(SiteCode) %>% 
  mutate(Area = median(Area,na.rm=T)) 

dat <- dat %>% 
  group_by(TRCA) %>% 
  mutate(Area_class = cut(`Usable Area`,
                          breaks = scales::log_breaks(5)(`Usable Area`)
  )) %>% 
  ungroup()


dat <- dat %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

levels(dat$`Stress Class`) <- c("Low","Medium","High")
dat$`Stress Class` <- as.character(dat$`Stress Class`)
dat$`Stress Class`[is.na(dat$`Stress Class`)] <- "North"
dat$`Stress Class` <- factor(dat$`Stress Class`,levels = c("Low","Medium","High","North"))

dat <- dat %>% 
  group_by(SiteCode) %>% 
  mutate(`Habitat Stability Class` = median(`Habitat Stability Score`,na.rm=T)) %>% 
  ungroup() %>% 
  mutate(`Habitat Stability Class` = cut(`Habitat Stability Class`,
                              breaks = 2)) 

levels(dat$`Habitat Stability Class`) <- c("Low","High")
dat$`Habitat Stability Class` <- as.character(dat$`Habitat Stability Class`)
dat$`Habitat Stability Class`[is.na(dat$`Habitat Stability Class`)] <- "Missing"
dat$`Habitat Stability Class` <- factor(dat$`Habitat Stability Class`,levels = c("Low","High","Missing"))

dat$`UCA Class` <- cut(dat$`Upstream Catchment Area`,c(0,20,80,Inf))
levels(dat$`UCA Class`) <- c("Small","Intermediate","Large")

# Pivot data long ---------------------------------------------------------

dat_long <- dat %>% 
  pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
               names_to = "Endpoint") 

# Separate locations ------------------------------------------------------

locs <- dat %>% 
  select(all_of(c(ID_cols[!ID_cols %in% c("ID","SampleDate")],"Sample Years","Years of Data","Sample Months","Sample Year-Months","Longitude","Latitude"))) %>% 
  distinct() %>% 
  sf::st_as_sf(coords=c("Longitude","Latitude"),crs="WGS84")

if (F) mapview::mapview(locs,zcol="OrgNameList")


# Save Data ---------------------------------------------------------------

saveRDS(
  list(
    wide = dat,
    long = dat_long
  ),
  "data/1_formatted_data.RDS"
)

# Calculate Fish Endpoints ------------------------------------------------


lookup_tbl<-read.csv(file.path("data","fish_lookup.csv")) %>% 
  mutate(name.to.use=trimws(name.to.use))
colnames(lookup_tbl)[colnames(lookup_tbl)=="OMNR.Code"]<-"SpeciesCode"

SATI_df<-read.csv(file.path("data","lkpSpeciesHabitatRequirement.csv")) %>% 
  select(SpeciesCode=FWISFishCode,OPT_TEMP) %>% 
  distinct() %>% 
  filter(!is.na(OPT_TEMP))

lookup_tbl<-left_join(lookup_tbl,SATI_df,by="SpeciesCode")

raw_tbl<-read.csv(file.path("data","tblFishSummaryOfTotalCatches.csv")) %>% 
  filter(SiteCode %in% locs$SiteCode,
         year(as.Date(SampleDate))>1995)

raw_tbl<-raw_tbl %>% 
  filter( 
    OSAPSE==1, #the sample event was associated with at least one OSAP project and the sample event itself used site boundaries that were defined as per OSAP
    !is.na(TotalWeightPer100m2) # This makes sure only valid single taxa are included (i.e., Cyprid bulk samples will be excluded)
  )

# Setup Sample Event table ------------------------------------------------

SampleEventID<-raw_tbl[!is.na(raw_tbl$UsableArea) & raw_tbl$UsableArea>0,
                       c("SampleEventID",
                         "StreamName",
                         "StreamCode",
                         "SiteCode",
                         "Latitude",
                         "Longitude",
                         "SampleEventType",
                         "SampleDate")]

SampleEventID<-distinct(SampleEventID)

SampleEventID<-SampleEventID[SampleEventID$SampleEventType=="(Electro) Fish Sampling",]
SampleEventID<-SampleEventID[!is.na(SampleEventID$Latitude) & !is.na(SampleEventID$Longitude),]

# Setup Fish table --------------------------------------------------------

sub_tbl<-raw_tbl[raw_tbl$SampleEventID %in% SampleEventID$SampleEventID,
                 c("SampleEventID","SampleDate","SpeciesCode","CommonName","ScientificName","NumberOfFishPer100m2","TotalWeightPer100m2")] 

# sub_tbl$EstimatedBiomass<-sub_tbl$BulkWeight/sub_tbl$UsableArea*100
# sub_tbl$NumberOfFish<-sub_tbl$NumberOfFish/sub_tbl$UsableArea*100
sub_tbl$EstimatedBiomass<-sub_tbl$TotalWeightPer100m2
sub_tbl$NumberOfFish<-sub_tbl$NumberOfFishPer100m2

sub_tbl<-sub_tbl[sub_tbl$NumberOfFish>0,] %>% 
  group_by(SampleEventID, SampleDate, SpeciesCode,CommonName,ScientificName) %>% 
  summarise(
    NumberOfFish=sum(NumberOfFish,na.rm=T),
    EstimatedBiomass=sum(EstimatedBiomass,na.rm=T)
  ) %>% 
  ungroup() %>% 
  na.omit()

atr_tbl<-left_join(sub_tbl,lookup_tbl,by="SpeciesCode") 



# Taxa Tables -------------------------------------------------------------

taxa_tbl1 <- atr_tbl %>% 
  group_by(SampleEventID,SampleDate) %>% 
  mutate(Comm_Biomass=sum(EstimatedBiomass,na.rm=T),
         Comm_Abundance=sum(NumberOfFish,na.rm=T),
         Perc_Biomass=EstimatedBiomass/Comm_Biomass,
         Perc_Abundance=NumberOfFish/Comm_Abundance
  ) %>% 
  ungroup() %>% 
  select(SampleEventID,SampleDate,SpeciesCode,Comm_Biomass,Comm_Abundance,Perc_Biomass,Perc_Abundance) %>% 
  mutate(SpeciesCode=factor(SpeciesCode)) %>% 
  group_by(SampleEventID,SampleDate) %>% 
  complete(SpeciesCode,fill = list(Comm_Biomass=0,Comm_Abundance=0,Perc_Biomass=0,Perc_Abundance=0)) %>% 
  ungroup() %>% 
  mutate(Comm_Biomass=ifelse(Perc_Abundance>0 & Comm_Biomass==0,NA_real_,Comm_Biomass)) %>% 
  mutate(Perc_Biomass=ifelse(Perc_Abundance>0 & Perc_Biomass==0,NA_real_,Perc_Biomass))


# No Catch Table ----------------------------------------------------------

no_catch_tbl<-SampleEventID %>% 
  as_tibble() %>% 
  select(SampleEventID,SampleDate) %>% 
  left_join(taxa_tbl1,by = c("SampleEventID", "SampleDate")) %>% 
  filter(is.na(Comm_Biomass)) %>% 
  mutate_at(vars(Comm_Biomass:Perc_Abundance),~0) %>% 
  mutate(SpeciesCode=atr_tbl$SpeciesCode[[1]]) %>% 
  mutate(SpeciesCode=factor(SpeciesCode,levels=unique(taxa_tbl1$SpeciesCode))) %>% 
  group_by(SampleEventID,SampleDate) %>% 
  #complete(SpeciesCode,fill = list(Comm_Biomass=0,Comm_Abundance=0,Perc_Biomass=0,Perc_Abundance=0)) %>% 
  ungroup() 

lookup_tbl2<-read_csv(file.path("data","fish_lookup.csv")) %>% 
  mutate(`name to use`=trimws(`name to use`))
colnames(lookup_tbl2)[colnames(lookup_tbl2)=="OMNR Code"]<-"SpeciesCode"

taxa_tbl<-taxa_tbl1 %>% 
  bind_rows(no_catch_tbl) %>% 
  mutate(SpeciesCode=as.character(SpeciesCode)) %>% 
  mutate(SpeciesCode=as.integer(SpeciesCode)) %>% 
  left_join(lookup_tbl2,by="SpeciesCode") %>% 
  group_by(SampleEventID,SampleDate,`name to use`) %>% 
  reframe(across(where(is.numeric),~sum(.,na.rm=T)),
          across(where(is.character),~head(.[!is.na(.)],1))
  ) %>% 
  ungroup() %>% 
  select(any_of(colnames(lookup_tbl2)),everything())

# Calculate endpoints -----------------------------------------------------

atr_tbl<-left_join(sub_tbl,lookup_tbl,by="SpeciesCode") 

ep_tbl<-calc_fish_ep(atr_tbl)

# No Catch Table ----------------------------------------------------------

no_catch_tbl<-SampleEventID %>% 
  select(SampleEventID,SampleDate) %>% 
  left_join(ep_tbl,by = c("SampleEventID", "SampleDate")) %>% 
  filter(is.na(Comm_Richness)) %>% 
  mutate_at(vars(Comm_Biomass:Rep_Lith_Rich),~0)

ep_tbl[sapply(ep_tbl,is.nan)]<-NA

ep_tbl<-ep_tbl %>% 
  bind_rows(no_catch_tbl)


# Endpoint Index ----------------------------------------------------------
ep_tbl_lng<-ep_tbl %>% 
  gather("Endpoint","Value",-SampleEventID,-SampleDate) %>% 
  mutate(Value=ifelse(is.nan(Value),NA,Value)) %>% 
  mutate(Endpoint_group=case_when(
    grepl("^Therm_",Endpoint) ~ "Thermal Regime",
    grepl("^Troph_",Endpoint) ~ "Trophic Guild",
    grepl("^Tol_",Endpoint) ~ "Tolerance Group",
    grepl("^Spec_",Endpoint) ~ "Individual Species",
    grepl("^Env_",Endpoint) ~ "Environment Preference",
    grepl("^Rep_",Endpoint) ~ "Reproductive Guild",
    grepl("^Comm_",Endpoint) ~ "General Community"
  )) %>%
  mutate(Endpoint_subgroup=str_split(Endpoint,"_",simplify = T)[,2]) %>% 
  mutate(Endpoint_subgroup=case_when(
    Endpoint_subgroup == "SATI" ~ "SATI",
    Endpoint_subgroup == "cold" ~ "Cold water",
    Endpoint_subgroup == "cool" ~ "Cool water",
    Endpoint_subgroup == "coolcold" ~ "Cold-Cool water",
    Endpoint_subgroup == "warm" ~ "Warm water",
    Endpoint_subgroup == "invert" ~ "Invertivore",
    Endpoint_subgroup == "detrit" ~ "Detritivore",
    Endpoint_subgroup == "carn" ~ "Carniivore",
    Endpoint_subgroup == "herb" ~ "Herbiivore",
    Endpoint_subgroup == "plank" ~ "Planktivore",
    Endpoint_subgroup == "tol" ~ "Tolerant",
    Endpoint_subgroup == "intol" ~ "Intolerant",
    Endpoint_subgroup == "RBD" ~ "Rainbow Darter",
    Endpoint_subgroup == "BKT" ~ "Brook Trout",
    Endpoint_subgroup == "GDF" ~ "Goldfish",
    Endpoint_subgroup == "RDG" ~ "Round Goby",
    Endpoint_subgroup == "CKC" ~ "Creek Chub",
    Endpoint_subgroup == "CMC" ~ "Common Carp",
    Endpoint_subgroup == "SKP" ~ "Mottled and Slimy Sculpin",
    Endpoint_subgroup == "LMP" ~ "Lampreys",
    Endpoint_subgroup == "DRT" ~ "Sensitive Darters",
    Endpoint_subgroup == "CEB" ~ "Centrarchids and Basses",
    Endpoint_subgroup == "CTF" ~ "Catfishes",
    Endpoint_subgroup == "SAL" ~ "Salmonids",
    Endpoint_subgroup == "benthic" ~ "Benthic",
    Endpoint_subgroup == "pelagic" ~ "Pelagic",
    Endpoint_subgroup == "benthopelagic" ~ "Benthopelagic",
    Endpoint_subgroup == "guard" ~ "Guarders",
    Endpoint_subgroup == "BrHid" ~ "Brood Hiders/Nest Spawners",
    Endpoint_subgroup == "OSS" ~ "Open Substratum/Substratum Choosers",
    Endpoint_subgroup == "Lith" ~ "Lithophilic Spawners",
    Endpoint=="Comm_Biomass" ~ "Total Biomass",
    Endpoint=="Comm_Abundance" ~ "Total Abundance",
    Endpoint=="Comm_Richness" ~ "Total Richness"
  )) %>% 
  mutate(Endpoint_type=case_when(
    grepl("_Bioperc$",Endpoint) ~ "Biomass (%)",
    grepl("_Abuperc$",Endpoint) ~ "Abundance (%)",
    grepl("_Bio$",Endpoint) ~ "Total Biomass",
    grepl("_Abu$",Endpoint) ~ "Total Abundance",
    grepl("_Rich$",Endpoint) ~ "Richness",
    Endpoint=="Comm_Biomass" ~ "Total Biomass",
    Endpoint=="Comm_Abundance" ~ "Total Abundance",
    Endpoint=="Comm_Richness" ~ "Richness"
  )) %>% 
  filter(!is.na(Value))


# Write results to excel --------------------------------------------------

ep_tbl_out<-ep_tbl %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())

if (F){
  ep_tbl_out %>% 
    mutate(Year = year(as.Date(SampleDate)),
           month = month(as.Date(SampleDate))) %>% 
    group_by(SiteCode,Year) %>% 
    summarise(
      n_samples = n(),
      months = paste0(month,collapse = ", ")
    ) %>% 
    arrange(desc(n_samples)) %>% 
    View()
}

ep_tbl_lng_out<-ep_tbl_lng %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())

param_tbl_out<-ep_tbl_lng %>% 
  select(Endpoint_group,Endpoint_subgroup,Endpoint_type,Endpoint) %>% 
  distinct()

taxa_tbl_out<-taxa_tbl %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())


wb = createWorkbook()

addWorksheet(wb, "Endpoints")
addWorksheet(wb, "Results - Wide")
addWorksheet(wb, "Results - Long")
addWorksheet(wb, "Taxa Table")

writeData(wb,param_tbl_out, sheet="Endpoints", rowNames =FALSE)
writeData(wb,ep_tbl_out, sheet="Results - Wide", rowNames =FALSE)
writeData(wb,ep_tbl_lng_out, sheet="Results - Long", rowNames =FALSE)
writeData(wb,taxa_tbl_out, sheet="Taxa Table", rowNames =FALSE)

saveWorkbook(wb, file.path("data","Fish_Endpoints.xlsx"),overwrite = T)

final_fish_out <- ep_tbl_lng_out %>% 
  mutate(Year = year(as.Date(SampleDate))) %>% 
  select(SiteCode,Year,Endpoint,Endpoint_group,Endpoint_subgroup,Endpoint_type,Value) %>% 
  group_by(SiteCode,Year,Endpoint,Endpoint_group,Endpoint_subgroup,Endpoint_type) %>% 
  summarise(Value = median(Value,na.rm = T), .groups="drop") 

write_csv(final_fish_out,"data/fish_endpoints.csv")

4 Fit Models and Format Results

source("00_source.R")
library(glmmTMB)
library(DHARMa)
library(broom.mixed)
library(parameters)

dat_in <- readRDS("data/1_formatted_data.RDS")$long 
predictor_summaries<-read_csv("data/predictor_summaries.csv")

if (F) {
  ggplot(dat_in,aes(x=value)) +
    geom_histogram() +
    facet_wrap(~Endpoint,scales="free")
  
  ggplot(dat_in,aes(x=value)) +
    geom_histogram() +
    scale_x_continuous(trans="log1p")+
    facet_wrap(~Endpoint,scales="free")
}

# Setup Model param -------------------------------------------------------

dat_in <- dat_in %>% 
  filter(`Stress Class` != "North") %>% # Exclude sites without land cover data
  filter(Endpoint != c("Habitat Stability Score")) %>%  # Has a problematic distribution
  filter(Endpoint != c("Sediment Transport Potential")) %>%  # Has a problematic distribution
  group_by(TRCA) %>% 
  mutate(
    Year_n2 = Year,
  ) %>% 
  mutate(
    Year_n = (max(Year)-min(Year))-(max(Year)-Year)+1,
  ) %>% 
  ungroup() %>% 
  mutate(
    `Strahler Order` = as.factor(`Strahler Order`)
  )

ep_tbl <- enframe(met_cols) %>% 
  unnest(value) %>% 
  mutate(Endpoint = names(value)) %>% 
  unnest(value)

mod_data_fin <- dat_in %>% 
  as_tibble() %>% 
  group_by(across(any_of(c("TRCA","Family","Endpoint")))) %>% 
  nest() %>% 
  ungroup() %>% 
  left_join(enframe(met_cols) %>% 
              unnest(value) %>% 
              mutate(Endpoint = names(value)) %>% 
              unnest(value) %>% 
              select(-value,Endpoint_group=name,Endpoint),
            by = "Endpoint" ) %>% 
  mutate(Endpoint = factor(Endpoint,levels = ep_tbl$Endpoint),
         Endpoint_group = factor(Endpoint_group,levels = unique(ep_tbl$name))) %>% 
  mutate(model_params = case_when(
    Endpoint %in% c("D15",
                    "D50",
                    "D85",
                    "D50Max"
    ) ~ list(
      list(
        family = lognormal(link = "log"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Sediment Sorting Index"
    ) ~ list(
      list(
        family = ziGamma(link = "log"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Width:Depth Ratio",
                    "Width - Mean",
                    "Depth - Mean") ~ list(
                      list(
                        family = lognormal(link = "log"),
                        ziformula = ~ 0,
                        dispformula = ~ 1
                      )
                    ),
    Endpoint %in% c("Sediment Transport Potential") ~ list(
      list(
        family = Gamma(link = "sqrt"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Usable Area") ~ list(
      list(
        family = Gamma(link = "log"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    grepl(" Cover| Percent",Endpoint) | Endpoint %in% c("Habitat Homogeneity",
                                                        "Habitat Stability Score")  ~ list(
                                                          list(
                                                            family = beta_family(link = "logit"),
                                                            ziformula = ~ 1,
                                                            dispformula = ~ 1
                                                          )
                                                        ),
    Endpoint %in% c("Bank Stability - Mean")  ~ list(
      list(
        family = beta_family(link = "probit"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("") ~ list(
      list(
        family = gaussian(link = "identity"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    )
  )) %>% 
  mutate(model_params = map2(TRCA,model_params,
                             ~case_when(.x==T ~ c(.y,
                                                  list(
                                                    control = glmmTMBControl(parallel = parallel::detectCores()),
                                                    formula = as.formula("value ~ `Upstream catchment developed (%)` + `Upstream catchment agricultural (%)` + log1p(`Reach Slope`) + Year_n + `UCA Class` + `Habitat Stability Class` + `Stress Class` + Year_n:`UCA Class` + Year_n:`Habitat Stability Class` + Year_n:`Stress Class` + Year_n:`UCA Class`:`Stress Class` + Year_n:`UCA Class`:`Habitat Stability Class` + (Year_n|`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode)")
                                                  )),
                                        .x==F ~ c(.y,
                                                  list(
                                                    control = glmmTMBControl(parallel = parallel::detectCores()),
                                                    formula = as.formula("value ~ `Upstream catchment developed (%)` + `Upstream catchment agricultural (%)` + log1p(`Reach Slope`) + Year_n + `UCA Class` + `Habitat Stability Class` + `Stress Class` + Year_n:`UCA Class` + Year_n:`Habitat Stability Class` + Year_n:`Stress Class` + Year_n:`UCA Class`:`Stress Class` + Year_n:`UCA Class`:`Habitat Stability Class` + (Year_n|`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode)")
                                                  ))
                             )
                             
  )) %>% 
  mutate(model_params = map2(TRCA,model_params,
                             ~{
                               if (.y$dispformula[[2]]!=0 & .x==T){
                                 .y$dispformula <- as.formula(" ~ 1")
                               }
                               if (.y$dispformula[[2]]!=0 & .x==F){
                                 .y$dispformula <- as.formula(" ~ 1")
                               }
                               if (.y$ziformula[[2]]!=0 & .x==T){
                                 .y$ziformula <- as.formula(" ~ `Upstream catchment developed (%)` + `Upstream catchment agricultural (%)` + log1p(`Reach Slope`) + Year_n + `UCA Class` + `Habitat Stability Class` + `Stress Class` + Year_n:`UCA Class` + Year_n:`Habitat Stability Class` + Year_n:`Stress Class` + Year_n:`UCA Class`:`Stress Class` + Year_n:`UCA Class`:`Habitat Stability Class` + (Year_n|`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode)")
                               }
                               if (.y$ziformula[[2]]!=0 & .x==F){
                                 .y$ziformula <- as.formula(" ~ `Upstream catchment developed (%)` + `Upstream catchment agricultural (%)` + log1p(`Reach Slope`) + Year_n + `UCA Class` + `Habitat Stability Class` + `Stress Class` + Year_n:`UCA Class` + Year_n:`Habitat Stability Class` + Year_n:`Stress Class` + Year_n:`UCA Class`:`Stress Class` + Year_n:`UCA Class`:`Habitat Stability Class` + (Year_n|`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode)")
                               }
                               return(.y)
                             }
  )) %>% 
  mutate(model_params = map2(model_params,data, 
                             ~c(.x,
                                list(
                                  data= .y %>% 
                                    filter(!is.na(value)) %>% 
                                    select(
                                      any_of(
                                        c("WSHED","Subshed","SiteCode","Area_class",
                                          "Sample_Season","Year_n","Year_n2","value",
                                          "year_class","Year",
                                          "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                                          "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                                          "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class")
                                      )
                                    ) 
                                )
                             ))) %>% 
  mutate(Endpoint = factor(Endpoint, levels = unlist(sapply(met_cols,names),use.names=F)))


mod_data_fin <- mod_data_fin %>% 
  mutate(p0 = map_dbl(data,~sum(.x$value==0,na.rm=T)/length(.x$value==0)))

mod_data_fin <- mod_data_fin %>% 
  filter(p0 < 0.75)

# Run Models --------------------------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(
    model = map(model_params, ~do.call(glmmTMB::glmmTMB,.x))
  ) 

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map_lgl(model,
                         ~any(
                           is.na(
                             predict(.x,
                                     type="link",
                                     re.form = NULL,
                                     se.fit = TRUE,
                                     cov.fit = F,
                                     do.bias.correct = F,
                                     allow.new.levels = T)$se.fit
                           ) |
                             any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                         ) 
    )
  )

# Exclude outliers --------------------------------------------------------
mod_data_fin <- mod_data_fin %>% 
  mutate(
    model_params = case_when(
      pred_error == T ~ map2(model_params,model, 
                             ~{
                               .x$data <- .x$data[!DHARMa::outliers(DHARMa::simulateResiduals(.y, use.u = T,refit = F,plot = F,n=9999),return = "logical",lowerQuantile = 0.025, upperQuantile = 0.975),]
                               return(.x)
                             }
      ),
      T ~ model_params
    )
  ) 

# Refit Models --------------------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    model = case_when(
      pred_error == T ~ map(model_params, ~do.call(glmmTMB::glmmTMB,.x)), #refit models with default links
      T ~ model
    )
  )

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map_lgl(model,
                         ~any(
                           is.na(
                             predict(.x,
                                     type="link",
                                     re.form = NULL,
                                     se.fit = TRUE,
                                     cov.fit = F,
                                     do.bias.correct = F,
                                     allow.new.levels = T)$se.fit
                           ) |
                             any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                         ) 
    )
  )



print(mod_data_fin,n=Inf)

# adjust link functions for bad fits
mod_data_fin <- mod_data_fin %>% 
  mutate(model = pmap(
    list(.x=model,
         .y=pred_error,
         .z=model_params),
    function(.x,.y,.z){
      if (!.y) return(.x)
      #if (.x$call$family$family != "beta") return(.x)
      if (.x$call$family$family == "beta") {
        .z$family <- beta_family(link = "probit")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      if (.x$call$family$family == "Gamma") {
        .z$family <- ziGamma(link = "sqrt")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      if (.x$call$family$family == "lognormal") {
        .z$family <- lognormal(link = "sqrt")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      return(.x)
    })
  )

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map2_lgl(model,
                          pred_error,
                          ~{
                            if (!.y) return(.y)
                            any(
                              is.na(
                                predict(.x,
                                        type="link",
                                        re.form = NULL,
                                        se.fit = TRUE,
                                        cov.fit = F,
                                        do.bias.correct = F,
                                        allow.new.levels = T)$se.fit
                              )
                            ) |
                              any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                          })
  )

print(mod_data_fin,n=Inf)



mod_data_fin <- mod_data_fin %>% 
  filter(!pred_error)

# Predict observations ----------------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(Pred1 = map2(model,model_params, function (i,nm) {
    
    year_vals <- nm$data %>% 
      select(Year_n,Year_n2,year_class) %>% 
      distinct()  
    
    pred_frame <- nm$data %>% 
      select(any_of(c("WSHED","Subshed","Area_class","SiteCode",
                      "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                      "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                      "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class"))) %>% 
      distinct() %>% 
      mutate(Sample_Season="Summer",#list(factor(levels(nm$data$Sample_Season),levels = levels(nm$data$Sample_Season))),
             Year_n = list(seq(from=min(year_vals$Year_n),to=max(year_vals$Year_n)))
      ) %>%
      unnest(Year_n) %>%
      unnest(Sample_Season) %>% 
      left_join(
        year_vals
      )
    
    pred_frame<-pred_frame[complete.cases(pred_frame),]
    
    a1 <- pred_frame %>% 
      bind_cols(
        data.frame(
          pred=predict(i,
                       newdata = pred_frame,
                       type="link",
                       re.form = NULL,
                       se.fit = TRUE,
                       cov.fit = F,
                       do.bias.correct = F,
                       allow.new.levels = T)
        )
      ) %>%
      mutate(pred.fit = predict(i,newdata = pred_frame,type="response")) %>% 
      mutate(pred.fit = i$modelInfo$family$linkfun(pred.fit)) %>% 
      mutate(
        pred_low = pred.fit - (qnorm(1 - 0.1 / 2) * pred.se.fit),
        pred_high = pred.fit + (qnorm(1 - 0.1 / 2) * pred.se.fit),
        pred_low_trend = pred.fit - (qnorm(1 - 0.5 / 2) * pred.se.fit),
        pred_high_trend = pred.fit + (qnorm(1 - 0.5 / 2) * pred.se.fit),
      ) %>% 
      left_join(
        nm$data %>% select(any_of(c("WSHED","Subshed","SiteCode","Area_class",
                                    "Sample_Season","Year_n","Year_n2",
                                    "value","year_class",
                                    "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                                    "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                                    "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class")))
      ) %>% 
      mutate(
        across(
          c(pred.fit,pred_low,pred_high,pred_low_trend,pred_high_trend),
          ~i$modelInfo$family$linkinv(.x)
        )
      ) 
    
    return(a1)
  }))

mod_data_fin <- mod_data_fin %>% 
  mutate(Pred1_error = map_lgl(Pred1,~any(is.na(.x$pred.se.fit))))

print(mod_data_fin,n=Inf)

# Extract random effects terms --------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(rand_ef = map(model,~get_terms2(.x))) 


mod_data_fin <- mod_data_fin %>% 
  mutate(rand_ef_error = map_lgl(rand_ef,~any(is.na(.x$sd))))

print(mod_data_fin,n=Inf)

# Get simulated residuals -------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    sim_resid_refitF = map(model, ~ DHARMa::simulateResiduals(.x, use.u = T,refit = F,plot = F,n=9999)),
  )

# Format final data tables ------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    Pred2 = map(Pred1,
                ~ .x %>% 
                  group_by(WSHED, Subshed, Sample_Season, SiteCode, year_class,
                           `Upstream Catchment Area`,`Baseflow Index`,`Reach Slope`,`Strahler Order`,
                           `Upstream catchment developed (%)`,`Upstream catchment agricultural (%)`,
                           `Upstream catchment stress index`,`Stress Class`,`Habitat Stability Class`,`UCA Class`) %>% 
                  summarise(
                    across(
                      c(pred.fit, pred.se.fit, pred_low, pred_high, pred_low_trend, pred_high_trend, value),
                      ~median(.x, na.rm=T)
                    ),
                    .groups = "drop"
                  ) %>% 
                  mutate(`UCA Class`=as.character(`UCA Class`))
    ),
    Pred2 = map2(Pred2,rand_ef,
                 function(x,y){
                   x %>% 
                     left_join(
                       y %>% 
                         filter(term == "Year_n") %>% 
                         filter(group_var== "`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode") %>% 
                         filter(component == "cond") %>% 
                         select(-group_var,-term,-sd,-`1`) %>% 
                         rename(lower_slope_sig = pred_low,
                                upper_slope_sig = pred_high,
                                lower_slope_trend = pred_low_trend,
                                upper_slope_trend = pred_high_trend,
                                mean_slope = value,
                         )
                     )
                 }),
    Pred2 = map(Pred2,
                ~.x %>% 
                  mutate(
                    linear_trend = case_when(
                      lower_slope_sig > 0 ~ "Increasing Significant",
                      upper_slope_sig < 0 ~ "Decreasing Significant",
                      lower_slope_trend > 0 ~ "Increasing Trend",
                      upper_slope_trend < 0 ~ "Decreasing Trend",
                      T ~ "No Change"
                    )
                  ) %>% 
                  mutate(
                    linear_trend = factor(linear_trend,
                                          levels = c(
                                            "Increasing Significant",
                                            "Increasing Trend",
                                            "No Change",
                                            "Decreasing Trend",
                                            "Decreasing Significant"
                                          ))
                  )
                
    )
  )


# UCA, Stability and stress class trends ------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(
    spatial_summ = map(data,~summary_fun(.x))
  ) %>% 
  mutate(summary_terms_temporal = map(rand_ef,
                                      ~{
                                        .x %>% 
                                          filter(term == "Year_n") %>% 
                                          filter(group_var!= "`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode") %>% 
                                          filter(component == "fixed") %>% 
                                          rename(lower_slope_sig = pred_low,
                                                 upper_slope_sig = pred_high,
                                                 lower_slope_trend = pred_low_trend,
                                                 upper_slope_trend = pred_high_trend,
                                                 mean_slope = value,
                                          ) %>% 
                                          mutate(
                                            linear_trend = case_when(
                                              lower_slope_sig > 0 ~ "Increasing Significant",
                                              upper_slope_sig < 0 ~ "Decreasing Significant",
                                              lower_slope_trend > 0 ~ "Increasing Trend",
                                              upper_slope_trend < 0 ~ "Decreasing Trend",
                                              T ~ "No Change"
                                            )
                                          ) %>% 
                                          mutate(
                                            linear_trend = factor(linear_trend,
                                                                  levels = c(
                                                                    "Increasing Significant",
                                                                    "Increasing Trend",
                                                                    "No Change",
                                                                    "Decreasing Trend",
                                                                    "Decreasing Significant"
                                                                  ))
                                          ) 
                                      })) %>% 
  mutate(summary_terms_spatial = map(rand_ef,
                                     ~{.x %>% 
                                         filter(term == "(Intercept)") %>% #labelled intercept but is predicted at max year
                                         filter(group_var != "`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode") %>% 
                                         filter(component == "fixed") 
                                     }))


saveRDS(mod_data_fin,"data/results_out.rds")

# Model Fit tables --------------------------------------------------------------

mod_tbl_list <- list()

mod_tbl_list$TRCA_mod_tbl <- modelsummary::modelsummary(
  setNames(
    mod_data_fin$model[mod_data_fin$TRCA],
    mod_data_fin$Endpoint[mod_data_fin$TRCA]
  ),
  vcov = sandwich::vcovHC,
  shape = effect:component + term + group ~ model + statistic,
  stars = TRUE,
  gof_omit = c("BIC|AIC|Std|Cond|Marg|ICC"),
  gof_function = function(model) tibble(`Pseudo R2`=scales::label_number(0.001)(cor(predict(model),model$frame$value,method="spearman")^2)) ,
  output = "gt",
  coef_rename = F,
  notes = "zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;",
  title = "TRCA Data Only"
) 

sig_list <- which(apply(mod_tbl_list$TRCA_mod_tbl$`_data`,1,function(x)grepl("\\*|\\+",x)),arr.ind = T)
sig_list <- split(sig_list,1:nrow(sig_list))

for (i in sig_list){
  mod_tbl_list$TRCA_mod_tbl <- mod_tbl_list$TRCA_mod_tbl %>% 
    gt::tab_style(
      style = gt::cell_fill(color = "lightblue"),
      locations = gt::cells_body(
        columns = i[[1]],
        rows = i[[2]]
      )
    )
}

if (F){
  mod_tbl_list$ALL_mod_tbl <- modelsummary::modelsummary(
    setNames(
      mod_data_fin$model[!mod_data_fin$TRCA],
      mod_data_fin$Endpoint[!mod_data_fin$TRCA]
    ),
    vcov = sandwich::vcovHC,
    shape = effect:component + term + group ~ model + statistic,
    stars = TRUE,
    gof_omit = c("BIC|AIC|Std|Cond|Marg|ICC"),
    gof_function = function(model) tibble(`Pseudo R2`=scales::label_number(0.001)(cor(predict(model),model$frame$value,method="spearman")^2)),
    output = "gt",
    coef_rename = F,
    notes = "zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;",
    title = "All Data"
  ) 
  
  sig_list <- which(apply(mod_tbl_list$ALL_mod_tbl$`_data`,1,function(x)grepl("\\*|\\+",x)),arr.ind = T)
  sig_list <- split(sig_list,1:nrow(sig_list))
  
  for (i in sig_list){
    mod_tbl_list$ALL_mod_tbl <- mod_tbl_list$ALL_mod_tbl %>% 
      gt::tab_style(
        style = gt::cell_fill(color = "lightblue"),
        locations = gt::cells_body(
          columns = i[[1]],
          rows = i[[2]]
        )
      )
  }
}


saveRDS(mod_tbl_list,"data/mod_tbl_list.rds")

# Correlation Figures -----------------------------------------------------

cor_list <- list(
  TRCA_spat = mod_data_fin %>% 
    filter(TRCA) %>% 
    select(Endpoint, data) %>% 
    unnest(data) %>% 
    pivot_wider(names_from = Endpoint,values_from = value) %>% 
    select(any_of(unlist(sapply(met_cols,names),use.names=F)),
           any_of(c("Upstream Catchment Area","Reach Slope","Upstream catchment developed (%)","Upstream catchment agricultural (%)","Upstream catchment stress index"))
    ) 
  ,
  All_spat = mod_data_fin %>% 
    filter(!TRCA) %>% 
    select(Endpoint, data) %>% 
    unnest(data) %>% 
    pivot_wider(names_from = Endpoint,values_from = value) %>% 
    select(any_of(unlist(sapply(met_cols,names),use.names=F)),
           any_of(c("Upstream Catchment Area","Reach Slope","Upstream catchment developed (%)","Upstream catchment agricultural (%)","Upstream catchment stress index"))
    )
)

cor_list <- map(cor_list,function(x){
  heatmap_fn(x,
             alpha=0.1,
             sig_est=0.6,
             method="spearman",
             scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                      mid = "white",
                                                      high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                      limits=c(-1,1))) 
})



trend_cor <- mod_data_fin %>% 
  filter(TRCA) %>% 
  select(Endpoint, Pred2) %>% 
  mutate(Pred_current = map(Pred2,
                            ~.x %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer"))) %>% 
  select(-Pred2) %>% 
  unnest(Pred_current) %>% 
  select(SiteCode,Endpoint,linear_trend ) %>% 
  mutate(linear_trend = case_when(
    linear_trend == "Increasing Significant" ~ 2,
    linear_trend == "Increasing Trend" ~ 1,
    linear_trend == "No Change" ~ 0,
    linear_trend == "Decreasing Trend" ~ -1,
    linear_trend == "Decreasing Significant" ~ -2,
  )) %>% 
  pivot_wider(values_from = linear_trend ,names_from = Endpoint) %>% 
  select(-SiteCode)

t1 <- cor(trend_cor,method = "spearman",use="pairwise.complete.obs")
keep <- (!apply(t1,2,function(x) all(is.na(x))))
trend_cor <- trend_cor[,keep]

trend_cor<-heatmap_fn(
  trend_cor,
  alpha=0.1,
  sig_est=0.6,
  method="spearman",
  use="pairwise.complete.obs",
  scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                           mid = "white",
                                           high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                           limits=c(-1,1))) 

cor_list$TRCA_temp <- trend_cor

trend_cor <- mod_data_fin %>% 
  filter(!TRCA) %>% 
  select(Endpoint, Pred2) %>% 
  mutate(Pred_current = map(Pred2,
                            ~.x %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer"))) %>% 
  select(-Pred2) %>% 
  unnest(Pred_current) %>% 
  select(SiteCode,Endpoint,linear_trend ) %>% 
  mutate(linear_trend = case_when(
    linear_trend == "Increasing Significant" ~ 2,
    linear_trend == "Increasing Trend" ~ 1,
    linear_trend == "No Change" ~ 0,
    linear_trend == "Decreasing Trend" ~ -1,
    linear_trend == "Decreasing Significant" ~ -2,
  )) %>% 
  pivot_wider(values_from = linear_trend ,names_from = Endpoint) %>% 
  select(-SiteCode)

t1 <- cor(trend_cor,method = "spearman",use="pairwise.complete.obs")
keep <- (!apply(t1,2,function(x) all(is.na(x))))
trend_cor <- trend_cor[,keep]

trend_cor<-heatmap_fn(
  trend_cor,
  alpha=0.1,
  sig_est=0.6,
  method="spearman",
  use="pairwise.complete.obs",
  scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                           mid = "white",
                                           high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                           limits=c(-1,1))) 

cor_list$All_temp <- trend_cor

saveRDS(cor_list,"data/cor_plots.rds")

5 Generate Figures

source("00_source.R")
library(mapview)
library(basemaps)
library(ggspatial)

cols_list <- tail(scales::col_darker(RColorBrewer::brewer.pal(4,"Set1"),5),3)
names(cols_list)<-c("Small","Intermediate","Large")

mod_data_fin <- readRDS("data/results_out.rds")
coords <- readRDS("data/1_formatted_data.RDS")$wide %>% 
  select(TRCA,WSHED,SiteCode,`Stress Class`,`UCA Class`,`Habitat Stability Class`) %>% 
  group_by(TRCA,WSHED,SiteCode,`Stress Class`,`UCA Class`,`Habitat Stability Class`) %>% 
  summarise(geometry = unique(st_centroid(geometry)),.groups="drop")

saveRDS(coords,"data/coords.rds")

coords<-coords %>%
  select(-TRCA,-WSHED,-`Stress Class`,-`UCA Class`,-`Habitat Stability Class`) %>% 
  distinct()

sf::sf_use_s2(F)

TRCA_subwats <- sf::read_sf("data/Subwatersheds_TRCA.geojson") 
TRCA_wats <- TRCA_subwats %>% 
  group_by(WSHED) %>% 
  summarise(geometry = sf::st_union(geometry))

unzip("data/AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario.zip",exdir=tempdir())
strm_lns <- sf::read_sf(
  file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w03_Lake_Ontario_West/w03_AEC_V3_Core_Rev0.gpkg"),
  "w03_Reach_SimpleGeometry_Core") %>% 
  sf::st_transform(sf::st_crs(TRCA_subwats)) 

strm_lns <- strm_lns[apply(sf::st_intersects(strm_lns,sf::st_buffer(TRCA_wats,-0.0001),sparse=F),1,any),]

strm_lns <- strm_lns %>% 
  filter(Network_Line_Type != "Shoreline Virtual Connector")

mod_data_fin <- mod_data_fin %>% 
  filter(TRCA)

plot_out <- mod_data_fin %>% 
  mutate(fam = map_chr(model_params,~.x$family$family)) %>% 
  mutate(trans = case_when(
    fam == "beta" ~ list(clogit),
    fam == "lognormal" ~ list(scales::transform_pseudo_log()),
    fam == "Gamma" ~ list(scales::transform_pseudo_log()),
  )) %>% 
  mutate(brks = case_when(
    fam == "beta" ~ list(scales::pretty_breaks),
    fam == "lognormal" ~ list(pseudo_log_breaks),
    fam == "Gamma" ~ list(pseudo_log_breaks),
  )) %>% 
  select(TRCA, Endpoint, Endpoint_group, fam, trans, brks, Pred1, Pred2, summary_terms_temporal, summary_terms_spatial) %>% 
  mutate(Pred_current = map(Pred2,
                            ~left_join(.x,coords,by = join_by(SiteCode)) %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer") %>% 
                              st_as_sf() %>% 
                              filter(is.na(as.numeric(as.character(WSHED)))) %>% 
                              mutate(`UCA Class` = factor(`UCA Class`,levels=c("Small","Intermediate","Large")))
  ))%>% 
  mutate(Pred_obs = map(Pred1,
                        ~.x %>% 
                          distinct() %>% 
                          filter(!is.na(value)) %>% 
                          mutate(Year_n = 2000+Year_n) %>% 
                          filter(is.na(as.numeric(as.character(WSHED)))) %>% 
                          mutate(`UCA Class` = factor(`UCA Class`,levels=c("Small","Intermediate","Large")))
  )) %>% 
  select(-Pred2,-Pred1) 


# Raw Data Plots ----------------------------------------------------------

plot_out$raw_plots <- pmap(
  list(.x=plot_out$Pred_obs,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) 
    .x %>% 
    mutate(WSHED2=WSHED) %>% 
    group_by(WSHED2) %>% 
    nest() %>% 
    ungroup() %>% 
    mutate(plt = map2(data,.y,
                      function(x,y) {
                        ggplot(
                          x,
                          aes(x = Year_n,
                              colour = SiteCode,
                              group = SiteCode
                          )
                        ) +
                          geom_linerange(aes(ymin=pred_low,ymax=pred_high),
                                         colour='grey',
                                         alpha=0.5,
                                         linewidth = 5) +
                          geom_point(aes(y=value)) +
                          scale_colour_discrete() +
                          theme_bw() +
                          ylab(paste0(y)) +
                          xlab("")+
                          facet_wrap(~SiteCode) +
                          labs(
                            title = paste0(x$WSHED[[1]]," | ",.y),
                            subtitle = "Points represent observed values, grey bars represent 90 percent confidence interval from GLMM model"
                          )+
                          theme(
                            axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
                            legend.position = "none"
                          )
                      })) %>% 
    pull(plt) %>% 
    map(~.x + scale_y_continuous(breaks = .a(), expand = c(0,0)))%>% 
    map(~.x + coord_transform(y = .z))
)

saveRDS(plot_out,"data/rawdata_plots.rds")

plot_out <- plot_out %>% 
  select(-raw_plots)

# Model fit plot ----------------------------------------------------------

plot_out$model_fit <- pmap(
  list(.x=plot_out$Pred_obs,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = pred.fit,
          y = value
      )
    ) +
      geom_point(colour='grey25') +
      geom_abline(slope = 1, intercept = 0, linewidth = 0.25,colour = "black")+
      theme_bw() +
      ylab("Observed") +
      xlab("Model Fitted") +
      scale_x_continuous(limits = range(.a(5)(c(.x$pred.fit,.x$value))),breaks = .a()) +
      scale_y_continuous(limits = range(.a(5)(c(.x$pred.fit,.x$value))),breaks = .a()) +
      coord_transform(
        x = .z,
        y = .z
      ) + 
      labs(
        title = paste0(.y),
        subtitle = paste0("Pseudo R2 = ", scales::label_number(0.001)(cor(.x$pred.fit,.x$value,method="spearman")^2))
      )
  }
)


# Model Residuals --------------------------------------------------------------

plot_out$fit1 <- map2(mod_data_fin$sim_resid_refitF,
                      mod_data_fin$Endpoint,
                      function(x,y){
                        ggplotify::as.ggplot(function() plot(x,title = paste("Model Fits -",y)))
                      }
)


# Temporal Plots ----------------------------------------------------------

plot_out$temporal_plot1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = mean_slope,
          colour = WSHED,
          alpha = case_when(
            grepl("Sig",linear_trend) ~ 1,
            grepl("Trend",linear_trend) ~ 0.92,
            T ~ 0.8
          )
      )
    ) +
      geom_pointrange(aes(ymin = lower_slope_sig,
                          ymax = upper_slope_sig,
      ),linewidth = 0.25) +
      geom_pointrange(aes(ymin = lower_slope_trend,
                          ymax = upper_slope_trend,
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = mean_slope),
                 shape=21,size=2.5,inherit.aes = F)+
      geom_hline(yintercept = 0) +
      scale_colour_discrete() +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "none"
      ) +
      labs(title = paste0(.y),
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) 
  }
)

plot_out$temporal_plot1_1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = mean_slope,
          alpha = case_when(
            grepl("Sig",linear_trend) ~ 1,
            grepl("Trend",linear_trend) ~ 0.92,
            T ~ 0.8
          )
      )
    ) +
      geom_linerange(aes(ymin = lower_slope_sig,
                         ymax = upper_slope_sig,
                         colour = `UCA Class`
      ),linewidth = 0.25) +
      geom_linerange(aes(ymin = lower_slope_trend,
                         ymax = upper_slope_trend,
                         colour = `UCA Class`
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = mean_slope,
                     shape = `Stress Class`,
                     fill = `UCA Class`,
                     alpha = case_when(
                       grepl("Sig",linear_trend) ~ 1,
                       grepl("Trend",linear_trend) ~ 0.92,
                       T ~ 0.8
                     )),
                 colour = "black",
                 size=2.5,inherit.aes = F)+
      geom_hline(yintercept = 0) +
      scale_shape_manual(values=c(21,22,23)) +
      scale_colour_manual(values = cols_list) +
      scale_fill_manual(values = cols_list) +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "top"
      ) +
      labs(title = paste0(.y),
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) +
      guides(alpha="none")
  }
)


# Temporal plots by station -----------------------------------------------

plot_out_stntemp <- plot_out %>% 
  select(TRCA,Endpoint,Endpoint_group,fam,trans,brks,Pred_current) %>% 
  unnest(Pred_current) %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(
    temporal_plot1_st = map(data,
                            ~{
                              ggplot(
                                .x 
                              ,
                                aes(x = Endpoint,
                                    y = mean_slope,
                                    alpha = case_when(
                                      grepl("Sig",linear_trend) ~ 1,
                                      grepl("Trend",linear_trend) ~ 0.92,
                                      T ~ 0.8
                                    )
                                )
                              ) +
                                geom_linerange(aes(ymin = lower_slope_sig,
                                                   ymax = upper_slope_sig,
                                                   colour = Endpoint_group
                                ),linewidth = 0.25) +
                                geom_linerange(aes(ymin = lower_slope_trend,
                                                   ymax = upper_slope_trend,
                                                   colour = Endpoint_group
                                ),linewidth = 1.5) +
                                geom_point(aes(x = Endpoint,
                                               y = mean_slope,
                                               fill = Endpoint_group,
                                               alpha = case_when(
                                                 grepl("Sig",linear_trend) ~ 1,
                                                 grepl("Trend",linear_trend) ~ 0.92,
                                                 T ~ 0.8
                                               )),
                                           colour = "black",
                                           shape=21,
                                           size=2.5,inherit.aes = F)+
                                geom_hline(yintercept = 0) +
                                facet_wrap(~WSHED+SiteCode, scales = "free",ncol=4) +
                                theme_bw() +
                                ylab("") +
                                xlab("")+
                                theme(
                                  axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
                                  legend.position = "none"
                                ) +
                                labs(#title = paste0(.y),
                                  subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
                                  caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
                                ) +
                                guides(alpha="none")
                              
                            })
  )

saveRDS(plot_out_stntemp,"data/plot_out_stntemp.rds")

# Temporal Plots - High Level ---------------------------------------------

plot_out$temporal_plot2 <- pmap(
  list(.x=plot_out$summary_terms_temporal,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    .x %>% 
      mutate(
        gp = case_when(
          group_var == "`1`" ~ "Overall",
          group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
          group_var == "`Stress Class`" ~ `Stress Class`,
          group_var == "`UCA Class`" ~ `UCA Class`,
        ),
        gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))
      ) %>% 
      filter(group_var == "`UCA Class`:`Stress Class`:`Habitat Stability Class`") %>% 
      mutate(`Stress Class` = paste0("Stress - ",`Stress Class`)) %>%
      mutate(`Habitat Stability Class` = paste0("Habitat Stability - ",`Habitat Stability Class`)) %>%
      mutate(`Stress Class` = factor(`Stress Class`,levels = paste0("Stress - ",c("High","Medium","Low")))) %>% 
      mutate(`Habitat Stability Class` = factor(`Habitat Stability Class`,levels = paste0("Habitat Stability - ",c("High","Low")))) %>% 
      mutate(gp=case_when(
        is.na(`Stress Class`) ~ `Habitat Stability Class`,
        T ~ `Stress Class`)) %>% 
      ggplot(
        aes(x = `UCA Class`,
            y = mean_slope,
            colour = `Habitat Stability Class`,
            group = `Habitat Stability Class`,
            alpha = case_when(
              grepl("Sig",linear_trend) ~ 1,
              grepl("Trend",linear_trend) ~ 0.92,
              T ~ 0.8
            )
        )
      ) +
      geom_pointrange(aes(ymin = lower_slope_sig,
                          ymax = upper_slope_sig,
      ),linewidth = 0.25,position=position_dodge(width=0.5)) +
      geom_pointrange(aes(ymin = lower_slope_trend,
                          ymax = upper_slope_trend,
      ),linewidth = 1.5,position=position_dodge(width=0.5)) +
      geom_point(aes(x = `UCA Class`,
                     y = mean_slope,
                     group = `Habitat Stability Class`),
                 shape=21,size=2.5,inherit.aes = F,position=position_dodge(width=0.5))+
      geom_hline(yintercept = 0) +
      scale_colour_discrete() +
      facet_grid(~`Stress Class`) +
      theme_bw() +
      ylab(.y) +
      theme(
        legend.position = "bottom"
      ) +
      labs(title = paste0(.y, " by Stress and Habitat Stability Class"),
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Coefficients with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) +
      guides(alpha="none")
    
  }
)


# Temporal Plots - Summary ------------------------------------------------

temp_summary_plot <- plot_out %>%
  select(TRCA,
         summary_terms_temporal,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var %in% c("`Habitat Stability Class`","`Stress Class`","`UCA Class`","`1`")) %>% 
  mutate(
    gp = case_when(
      group_var == "`1`" ~ "Overall",
      group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
      group_var == "`Stress Class`" ~ `Stress Class`,
      group_var == "`UCA Class`" ~ `UCA Class`,
    ),
    group_var = case_when(
      group_var == "`1`" ~ "Overall",
      group_var == "`Habitat Stability Class`" ~ "Habitat Stability Class",
      group_var == "`Stress Class`" ~ "Stress Class",
      group_var == "`UCA Class`" ~ "UCA Class",
      T ~ group_var
    )
  ) %>% 
  mutate(group_var = factor(group_var, levels = rev(c("Habitat Stability Class","Stress Class","UCA Class","Overall")))) %>% 
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         mutate(gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))) %>% 
                         #filter(!is.na(gp)) %>% 
                         ggplot(
                           aes(x = gp,
                               y = mean_slope,
                               colour = group_var,
                               alpha = case_when(
                                 grepl("Sig",linear_trend) ~ 1,
                                 grepl("Trend",linear_trend) ~ 0.92,
                                 T ~ 0.8
                               )
                           )
                         ) +
                         geom_pointrange(aes(ymin = lower_slope_sig,
                                             ymax = upper_slope_sig,
                         ),linewidth = 0.25) +
                         geom_pointrange(aes(ymin = lower_slope_trend,
                                             ymax = upper_slope_trend,
                         ),linewidth = 1.5) +
                         geom_point(aes(x = gp,
                                        y = mean_slope),
                                    shape=21,size=2.5,inherit.aes = F) +
                         geom_hline(yintercept = 0) +
                         scale_colour_discrete() +
                         facet_grid(~group_var,scales="free", drop=TRUE) +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "none"
                         ) +
                         labs(title = .y)+
                         guides(alpha="none")
                     }))

temp_summary_plot <- temp_summary_plot %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr"))) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid( 
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Trends",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap 0 are considered to increase or decrease significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(temp_summary_plot,"data/temp_summary_plot.rds")

# Temporal Maps -----------------------------------------------------------

plot_out$temporal_map1 <- map2(
  plot_out$Pred_current,
  plot_out$Endpoint,
  ~basemap_ggplot(sf::st_bbox(TRCA_subwats),map_service="carto",map_type="light",verbose=F) +
    theme_minimal()+
    theme(legend.position = "bottom",
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          legend.key.width = unit(dev.size()[1] / 15, "inches"),
          legend.key.height = unit(dev.size()[2] / 50, "inches"),
          plot.margin=unit(c(0,0,0,0),"cm")) +
    ylab("")+
    xlab("")+
    geom_sf(data=sf::st_transform(TRCA_subwats,crs="EPSG: 3857"),alpha=0.75)+
    geom_sf(data=sf::st_transform(TRCA_wats,crs="EPSG: 3857"),alpha=0,linewidth=1.25)+
    geom_sf(data=sf::st_transform(strm_lns,crs="EPSG: 3857"),alpha=0.25,colour="darkblue")+
    geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
            aes(colour=linear_trend),
            size=3,
            linewidth=0.05) +
    geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
            size=3,
            shape=21) +
    scale_colour_manual(values = c(
      `Increasing Significant`="darkolivegreen4",
      `Increasing Trend`="darkolivegreen2",
      `No Change`="black",
      `Decreasing Trend`="skyblue2",
      `Decreasing Significant`="skyblue4"
    )) +
    labs(
      title = .y,
      subtitle  = "Mean slope coefficients from GLMM models",
      colour = .y
    ) 
)



# Spatial Plots -----------------------------------------------------------


plot_out$spatial_plot1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = pred.fit,
          colour = WSHED
      )
    ) +
      geom_pointrange(aes(ymin = pred_low,
                          ymax = pred_high,
      ),linewidth = 0.25) +
      geom_pointrange(aes(ymin = pred_low_trend,
                          ymax = pred_high_trend,
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = pred.fit),shape=21,size=2.5,inherit.aes = F)+
      
      scale_colour_discrete() +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "none"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Areas with confidence intervals that don't overlap 0 are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)

plot_out$spatial_plot1_1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = pred.fit,
          shape = `Stress Class`,
          colour = `UCA Class`
      )
    ) +
      geom_linerange(aes(ymin = pred_low,
                         ymax = pred_high,
                         colour = `UCA Class`
      ),linewidth = 0.25) +
      geom_linerange(aes(ymin = pred_low_trend,
                         ymax = pred_high_trend,
                         colour = `UCA Class`
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = pred.fit,
                     shape = `Stress Class`,
                     fill = `UCA Class`),size=2.5,inherit.aes = F)+
      scale_shape_manual(values=c(21,22,23))+
      scale_colour_manual(values = cols_list) +
      scale_fill_manual(values = cols_list) +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "top"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Areas with confidence intervals that don't overlap 0 are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)


# Spatial Plot - High Level -----------------------------------------------


plot_out$spatial_plot2 <- pmap(
  list(.x=plot_out$summary_terms_spatial,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    .x %>% 
      mutate(
        gp = case_when(
          group_var == "`1`" ~ "Overall",
          group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
          group_var == "`Stress Class`" ~ `Stress Class`,
          group_var == "`UCA Class`" ~ `UCA Class`,
        ),
        gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))
      ) %>% 
      filter(group_var == "`UCA Class`:`Stress Class`:`Habitat Stability Class`") %>% 
      mutate(`Stress Class` = paste0("Stress - ",`Stress Class`)) %>%
      mutate(`Habitat Stability Class` = paste0("Habitat Stability - ",`Habitat Stability Class`)) %>%
      mutate(`Stress Class` = factor(`Stress Class`,levels = paste0("Stress - ",c("High","Medium","Low")))) %>% 
      mutate(`Habitat Stability Class` = factor(`Habitat Stability Class`,levels = paste0("Habitat Stability - ",c("High","Low")))) %>% 
      mutate(gp=case_when(
        is.na(`Stress Class`) ~ `Habitat Stability Class`,
        T ~ `Stress Class`)) %>% 
      ggplot(
        aes(x = `UCA Class`,
            y = value,
            colour = `Habitat Stability Class`,
            group = `Habitat Stability Class`
        )
      ) +
      geom_pointrange(aes(ymin = pred_low,
                          ymax = pred_high,
      ),linewidth = 0.25,position=position_dodge(width=0.5)) +
      geom_pointrange(aes(ymin = pred_low_trend,
                          ymax = pred_high_trend,
      ),linewidth = 1.5,position=position_dodge(width=0.5)) +
      geom_point(aes(x = `UCA Class`,
                     y = value,
                     group = `Habitat Stability Class`),
                 shape=21,size=2.5,inherit.aes = F,position=position_dodge(width=0.5))+
      
      scale_colour_discrete() +
      facet_grid(~`Stress Class`) +
      theme_bw() +
      ylab(.y) +
      theme(
        legend.position = "bottom"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Coefficients with confidence intervals that don't overlap are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)

# Spatial Plots - Summary ------------------------------------------------

spat_summary_plot <- plot_out %>%
  select(TRCA,
         summary_terms_spatial,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_spatial) %>% 
  filter(group_var %in% c("`Habitat Stability Class`","`Stress Class`","`UCA Class`","`1`")) %>% 
  mutate(
    gp = case_when(
      group_var == "`1`" ~ "Overall",
      group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
      group_var == "`Stress Class`" ~ `Stress Class`,
      group_var == "`UCA Class`" ~ `UCA Class`,
    ),
    group_var = case_when(
      group_var == "`1`" ~ "Overall",
      group_var == "`Habitat Stability Class`" ~ "Habitat Stability Class",
      group_var == "`Stress Class`" ~ "Stress Class",
      group_var == "`UCA Class`" ~ "UCA Class",
      T ~ group_var
    )
  ) %>% 
  mutate(group_var = factor(group_var, levels = rev(c("Habitat Stability Class","Stress Class","UCA Class","Overall")))) %>% 
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         mutate(gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))) %>% 
                         ggplot(
                           aes(x = gp,
                               y = value,
                               colour = group_var
                           )
                         ) +
                         geom_pointrange(aes(ymin = pred_low,
                                             ymax = pred_high,
                         ),linewidth = 0.25) +
                         geom_pointrange(aes(ymin = pred_low_trend,
                                             ymax = pred_high_trend,
                         ),linewidth = 1.5) +
                         geom_point(aes(x = gp,
                                        y = value),
                                    shape=21,size=2.5,inherit.aes = F) +
                         scale_colour_discrete() +
                         facet_grid(~group_var,scales="free", drop=TRUE) +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "none"
                         ) +
                         labs(title=.y)+
                         guides(alpha="none") +
                         scale_y_continuous(breaks = .x$brks[[1]]()) +
                         coord_transform(
                           y = .x$trans[[1]]
                         )
                     }))

spat_summary_plot <- spat_summary_plot %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr") )) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid(
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Means",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap are considered to differ significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(spat_summary_plot,"data/spat_summary_plot.rds")

# Spatial Maps -----------------------------------------------------------

plot_out$spatial_map1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    basemap_ggplot(sf::st_bbox(TRCA_subwats),map_service="carto",map_type="light",verbose=F) +
      theme_minimal()+
      theme(legend.position = "bottom",
            axis.text = element_blank(),
            axis.ticks = element_blank(),
            legend.key.width = unit(dev.size()[1] / 15, "inches"),
            legend.key.height = unit(dev.size()[2] / 50, "inches"),
            plot.margin=unit(c(0,0,0,0),"cm")) +
      ylab("")+
      xlab("")+
      geom_sf(data=sf::st_transform(TRCA_subwats,crs="EPSG: 3857"),alpha=0.75)+
      geom_sf(data=sf::st_transform(TRCA_wats,crs="EPSG: 3857"),alpha=0,linewidth=1.25)+
      geom_sf(data=sf::st_transform(strm_lns,crs="EPSG: 3857"),alpha=0.25,colour="darkblue")+
      geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
              aes(colour=pred.fit),
              size=3,
              linewidth=0.1) +
      scale_colour_continuous(palette = "Spectral",type="div",trans = .z,breaks=.a()) +
      geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
              shape=21,
              size=3) +
      labs(
        title = .y,
        subtitle = "Predicted 2024 mean values from GLMM models",
        colour = .y
      ) 
    
  }
)


saveRDS(plot_out,"data/results_plots.rds")



# Temporal Summary Tables ----------------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,Pred_current) %>% 
  unnest(Pred_current) %>% 
  mutate(`UCA Class`=factor(`UCA Class`,levels = c("Small","Intermediate","Large"))) %>% 
  mutate(`Stress Class`=factor(`Stress Class`,levels = c("Low","Medium","High"))) %>% 
  mutate(`Habitat Stability Class`=factor(`Habitat Stability Class`,levels = c("High","Low"))) %>% 
  arrange(`Stress Class`,`UCA Class`,`Habitat Stability Class`) %>% 
  mutate(value = paste0(mean_slope," (",lower_slope_sig," - ",upper_slope_sig,")")) %>% 
  select(Endpoint_group,Endpoint,WSHED,SiteCode,value,linear_trend,
         `Stress Class`,`UCA Class`,`Habitat Stability Class`
  )

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(WSHED,SiteCode,`Stress Class`,`UCA Class`,`Habitat Stability Class`) %>% 
  distinct() %>% 
  group_by(WSHED) %>% 
  group_split() %>% 
  map(~setNames(list(.x %>% select(-WSHED)),.x$WSHED[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-value,-Endpoint_group,-WSHED) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  mutate(SiteCode = as.character(SiteCode)) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

writexl::write_xlsx(shade_tbl,"data/changetbl1.xlsx")

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

num_tbl <- tbl_out %>% 
  select(-linear_trend,-Endpoint_group,-WSHED) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  mutate(SiteCode = as.character(SiteCode))

col_list <- list(`Increasing Significant` = "darkolivegreen4",
                 `Increasing Trend` = "darkolivegreen2",
                 `No Change` = "grey90",
                 `Decreasing Trend` = "skyblue2",
                 `Decreasing Significant` = "skyblue4"
)


col_list_leg <- list(`+ + = Increasing Significant` = "darkolivegreen4",
                     `+ = Increasing Trend` = "darkolivegreen2",
                     `0 = No Change` = "darkgrey",
                     `- = Decreasing Trend` = "skyblue2",
                     `- - = Decreasing Significant` = "skyblue4"
)

col_list_leg <- map2(names(col_list_leg),col_list,
                     ~paste0("<span style='color:",col2hex(.y),";'>",.x,"</span>"))

col_list_leg <- paste0(col_list_leg, collapse = " | ")


gt_list <- list()

for (v in names(site_grp)){
  
  shade_tbl_sub <- shade_tbl %>% 
    filter(SiteCode %in% site_grp[[v]]$SiteCode) %>%
    select(-`Stress Class`)
  shade_tbl_sym_sub <- shade_tbl_sym %>%
    filter(SiteCode %in% site_grp[[v]]$SiteCode) %>%
    select(-`Stress Class`)
  
  gt_tbl <- gt::gt(shade_tbl_sym_sub  ,
                   rowname_col = "SiteCode") %>% 
    gt::tab_header(
      title = v
    ) %>%
    gt::tab_stubhead(label = v) %>% 
    gt::tab_source_note(
      source_note = gt::md(col_list_leg)
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center", v_align = "middle"),
      locations = gt::cells_body(
        columns = everything(),
        rows = everything()
      )
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center",
                            v_align = "middle",
                            stretch = "extra-condensed",
                            size = "small"),
      locations = gt::cells_column_labels()
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center",
                            v_align = "middle",
                            stretch = "extra-condensed",
                            size = "small",
                            transform="capitalize"),
      locations = gt::cells_row_groups()
    )
  
  for (i in names(col_list)){
    for (ii in 1:ncol(shade_tbl_sub)){
      gt_tbl <- gt_tbl %>% 
        gt::tab_style(
          style = gt::cell_fill(color = col_list[[i]]),
          locations = gt::cells_body(
            columns = ii,
            rows = which(shade_tbl_sub[[ii]] == i)
          )
        )
    }
  }
  
  for (i in names(ep_grp)){
    gt_tbl <- gt_tbl %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
  }
  
  for (i in unique(site_grp[[v]]$`Stress Class`)){
    gt_tbl <- gt_tbl %>% gt::tab_row_group(label = paste0(i," Stess Class"),rows = as.character(site_grp[[v]]$SiteCode[site_grp[[v]]$`Stress Class`==i]))
  }
  
  gt_list[[v]] <- gt_tbl
}

gt_tbl <- gt::gt_group(
  .list = gt_list
)

saveRDS(gt_tbl,"data/summarytbl.rds")

# Trend Summary Table 1 ---------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var %in% c("`1`","`Habitat Stability Class`","`Stress Class`","`UCA Class`")) %>% 
  select(Endpoint_group,Endpoint,`Stress Class`,`UCA Class`,`Habitat Stability Class`,`1`,linear_trend) %>% 
  arrange(`Stress Class`,`UCA Class`,`Habitat Stability Class`,`1`) %>% 
  pivot_longer(c(`Stress Class`,`UCA Class`,`Habitat Stability Class`,`1`)) %>% 
  filter(!is.na(value))

levels(tbl_out$value)[levels(tbl_out$value)=="overall"]<-"Global Trend"

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(name,value) %>% 
  distinct() %>% 
  group_by(name) %>% 
  group_split() %>% 
  map(~setNames(list(.x$value),.x$name[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-Endpoint_group) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  arrange(name,value) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

site_grp <- shade_tbl %>% 
  mutate(value = row_number()) %>% 
  select(name,value) %>% 
  distinct() %>% 
  group_by(name) %>% 
  group_split() %>% 
  map(~setNames(list(.x$value),.x$name[[1]])) %>% 
  unlist(recursive = F)

site_grp<-rev(site_grp)
names(site_grp)[names(site_grp)=="1"] <- " "


shade_tbl <- shade_tbl %>% 
  select(-name) 

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

gt_tbl1 <- gt::gt(shade_tbl_sym,rowname_col = "value") %>% 
  gt::tab_stubhead(label = "Overall Trends") %>% 
  gt::tab_source_note(
    source_note = gt::md(col_list_leg)
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center", v_align = "middle"),
    locations = gt::cells_body(
      columns = everything(),
      rows = everything()
    )
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small"),
    locations = gt::cells_column_labels()
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small",
                          transform="capitalize"),
    locations = gt::cells_row_groups()
  )

for (i in names(col_list)){
  for (ii in 1:ncol(shade_tbl)){
    gt_tbl1 <- gt_tbl1 %>% 
      gt::tab_style(
        style = gt::cell_fill(color = col_list[[i]]),
        locations = gt::cells_body(
          columns = ii,
          rows = which(shade_tbl[[ii]] == i)
        )
      )
  }
}

for (i in names(site_grp)){
  gt_tbl1 <- gt_tbl1 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
}

for (i in names(ep_grp)){
  gt_tbl1 <- gt_tbl1 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
}

# Trend Summary Table 2 ---------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var == "`UCA Class`:`Stress Class`") %>% 
  select(Endpoint_group,Endpoint,`Stress Class`,`UCA Class`,linear_trend) %>% 
  filter(!is.na(`Stress Class`)) %>% 
  mutate(`UCA Class` = paste("UCA",`UCA Class`)) %>% 
  mutate(`UCA Class` = factor(`UCA Class`,levels=c("UCA Small","UCA Intermediate","UCA Large"))) %>% 
  arrange(`UCA Class`)

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(`Stress Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Stress Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Stress Class`[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-Endpoint_group) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  arrange(`Stress Class`,`UCA Class`) %>% 
  mutate(`UCA Class` = as.character(`UCA Class`)) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

site_grp <- shade_tbl %>% 
  mutate(`UCA Class` = row_number()) %>% 
  select(`Stress Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Stress Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Stress Class`[[1]])) %>% 
  unlist(recursive = F)

site_grp <- site_grp[c("Low","Medium","High")]
names(site_grp) <- c(c("Low Stress","Medium Stress","High Stress"))

writexl::write_xlsx(shade_tbl,"data/changetbl2.xlsx")

shade_tbl <- shade_tbl %>% 
  select(-`Stress Class`) 

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

gt_tbl2 <- gt::gt(shade_tbl_sym,rowname_col = "UCA Class") %>% 
  # gt::tab_header(
  #   title = "Temporal Changes across UCA Class"
  # ) %>% 
  gt::tab_stubhead(label = "UCA Class / Stress Class") %>% 
  gt::tab_source_note(
    source_note = gt::md(col_list_leg)
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center", v_align = "middle"),
    locations = gt::cells_body(
      columns = everything(),
      rows = everything()
    )
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small"),
    locations = gt::cells_column_labels()
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small",
                          transform="capitalize"),
    locations = gt::cells_row_groups()
  )

for (i in names(col_list)){
  for (ii in 1:ncol(shade_tbl)){
    gt_tbl2 <- gt_tbl2 %>% 
      gt::tab_style(
        style = gt::cell_fill(color = col_list[[i]]),
        locations = gt::cells_body(
          columns = ii,
          rows = which(shade_tbl[[ii]] == i)
        )
      )
  }
}

for (i in names(site_grp)){
  gt_tbl2 <- gt_tbl2 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
}

for (i in names(ep_grp)){
  gt_tbl2 <- gt_tbl2 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
}


# Trend Summary Table 3 ---------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var == "`UCA Class`:`Habitat Stability Class`") %>% 
  select(Endpoint_group,Endpoint,`Habitat Stability Class`,`UCA Class`,linear_trend) %>% 
  filter(!is.na(`Habitat Stability Class`)) %>% 
  mutate(`UCA Class` = paste("UCA",`UCA Class`)) %>% 
  mutate(`UCA Class` = factor(`UCA Class`,levels=c("UCA Small","UCA Intermediate","UCA Large"))) %>% 
  arrange(`UCA Class`)


ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(`Habitat Stability Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Habitat Stability Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Habitat Stability Class`[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-Endpoint_group) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  arrange(`Habitat Stability Class`,`UCA Class`) %>% 
  mutate(`UCA Class` = as.character(`UCA Class`)) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

site_grp <- shade_tbl %>% 
  mutate(`UCA Class` = row_number()) %>% 
  select(`Habitat Stability Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Habitat Stability Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Habitat Stability Class`[[1]])) %>% 
  unlist(recursive = F)

site_grp <- site_grp[c("Low","High")]
names(site_grp) <- c(c("Low Stability","High Stability"))

writexl::write_xlsx(shade_tbl,"data/changetbl3.xlsx")

shade_tbl <- shade_tbl %>% 
  select(-`Habitat Stability Class`)

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

gt_tbl3 <- gt::gt(shade_tbl_sym,rowname_col = "UCA Class") %>% 
  gt::tab_stubhead(label = "UCA Class / Habitat Stability Class") %>% 
  gt::tab_source_note(
    source_note = gt::md(col_list_leg)
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center", v_align = "middle"),
    locations = gt::cells_body(
      columns = everything(),
      rows = everything()
    )
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small"),
    locations = gt::cells_column_labels()
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small",
                          transform="capitalize"),
    locations = gt::cells_row_groups()
  )

for (i in names(col_list)){
  for (ii in 1:ncol(shade_tbl)){
    gt_tbl3 <- gt_tbl3 %>% 
      gt::tab_style(
        style = gt::cell_fill(color = col_list[[i]]),
        locations = gt::cells_body(
          columns = ii,
          rows = which(shade_tbl[[ii]] == i)
        )
      )
  }
}

for (i in names(site_grp)){
  gt_tbl3 <- gt_tbl3 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
}

for (i in names(ep_grp)){
  gt_tbl3 <- gt_tbl3 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
}

gt_tbl4 <- gt::gt_group(
  gt_tbl1,
  gt_tbl2,
  gt_tbl3
)

saveRDS(gt_tbl4,"data/summarytbl2.rds")

6 Biological Correlations

source("00_source.R")

mod_data_fin <- readRDS("data/results_out.rds")

fish_ep <- read_csv("data/fish_endpoints.csv") %>% 
  select(SiteCode,Year,Endpoint,Value) %>% 
  pivot_wider(names_from = Endpoint, values_from = Value) %>% 
  select(SiteCode,Year, any_of(sel_fish_ep))

cor_list <- list()
fish_plot_list <- list()

# TRCA Data ---------------------------------------------------------------

cmb_data_TRCA <- mod_data_fin %>% 
  filter(TRCA) %>% 
  select(Endpoint, data) %>% 
  unnest(data) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  left_join(fish_ep,by = join_by(SiteCode, Year))

hab_data_TRCA <- cmb_data_TRCA %>% 
  select(any_of(unlist(sapply(met_cols,names),use.names=F))#,
         # any_of(c("Upstream Catchment Area","Reach Slope",
         #          "Upstream catchment stress index","Upstream catchment developed (%)","Upstream catchment agricultural (%)"))
  )

fish_data_TRCA <- cmb_data_TRCA %>% 
  select(any_of(names(sel_fish_ep)))

cor_list$TRCA <- heatmap_fn(hab_data_TRCA,
                            fish_data_TRCA,
                            alpha=0.1,
                            sig_est=0.6,
                            method="spearman",
                            scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                                     mid = "white",
                                                                     high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                                     limits=c(-1,1))) 

fish_plot_list$TRCA <- cor(fish_data_TRCA,hab_data_TRCA,method="spearman",use="pairwise.complete.obs") %>% 
  as.data.frame() %>% 
  rownames_to_column("var1") %>% 
  pivot_longer(c(everything(),-var1),names_to = "var2") %>% 
  mutate(abs_cor = abs(value)) %>% 
  filter(abs_cor>0.3)%>% 
  arrange(desc(abs_cor)) %>% 
  mutate(plot = map2(
    var1,
    var2,
    ~tibble(
      SiteCode = cmb_data_TRCA$StreamCode,
      SampleDate = cmb_data_TRCA$SampleDate,
      var2 = cmb_data_TRCA[[.y]],
      var1 = fish_data_TRCA[[.x]]
    ) %>% 
      ggplot(aes(x=var2,y=var1))+
      geom_point()+
      geom_smooth(se=F)+
      ylab(.x)+
      xlab(.y)+
      theme_bw()
    
  )) %>% 
  mutate(plot = map2(plot,value,
                     ~.x + labs(
                       #title = "TRCA Data Only",
                       subtitle = paste0("Spearman Correlation = ",scales::label_number(0.001)(.y)))
  ))

# All Data ----------------------------------------------------------------

cmb_data_all <- mod_data_fin %>% 
  filter(!TRCA) %>% 
  select(Endpoint, data) %>% 
  unnest(data) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  left_join(fish_ep,by = join_by(SiteCode, Year))

hab_data_all <- cmb_data_all %>% 
  select(any_of(unlist(sapply(met_cols,names),use.names=F)),
         any_of(c("Upstream Catchment Area","Reach Slope",
                  "Upstream catchment stress index","Upstream catchment developed (%)","Upstream catchment agricultural (%)"))
  )

fish_data_all <- cmb_data_all %>% 
  select(any_of(names(sel_fish_ep)))

cor_list$all <- heatmap_fn(hab_data_all,
                            fish_data_all,
                            alpha=0.1,
                            sig_est=0.6,
                            method="spearman",
                            scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                                     mid = "white",
                                                                     high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                                     limits=c(-1,1))) 


fish_plot_list$all <- cor(fish_data_all,hab_data_all,method="spearman",use="pairwise.complete.obs") %>% 
  as.data.frame() %>% 
  rownames_to_column("var1") %>% 
  pivot_longer(c(everything(),-var1),names_to = "var2") %>% 
  mutate(abs_cor = abs(value)) %>% 
  filter(abs_cor>0.3)%>% 
  arrange(desc(abs_cor)) %>% 
  mutate(plot = map2(
    var1,
    var2,
    ~tibble(
      SiteCode = cmb_data_all$StreamCode,
      SampleDate = cmb_data_all$SampleDate,
      var2 = hab_data_all[[.y]],
      var1 = fish_data_all[[.x]]
    ) %>% 
      ggplot(aes(x=var2,y=var1))+
      geom_point()+
      geom_smooth(se=F)+
      ylab(.x)+
      xlab(.y)+
      theme_bw()
    
  )) %>% 
  mutate(plot = map2(plot,value,
                     ~.x + labs(
                       #title = "All Data",
                       subtitle = paste0("Spearman Correlation = ",scales::label_number(0.001)(.y)))
  ))

saveRDS(fish_plot_list,"data/biofish_plots.rds")
saveRDS(cor_list,"data/biocor_plots.rds")

Raw with Model Predictions Data

These figures show the raw data for each endpoint over time, grouped by watershed. The grey bars represent the 90 percent prediction interval for each observation.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Results - Model Summaries

Model Summaries

These tables show the model coefficients, estimated standard error, and significance for fixed effects terms. Random effects are summarized as variance and correlations within groups. Variance measures how much the random effects vary across different groups or levels in the hierarchy.Random effects models assume that observations within the same group are correlated. This correlation arises because these observations share common characteristics. Summary measures are provided below the model results.

TRCA Data Only

effect			Width - Mean		Depth - Mean		Width:Depth Ratio		Usable Area		D50		D85		D50Max		Sediment Sorting Index		Bank Stability - Mean		Habitat Homogeneity		Filamentous Algae Percent		Non-Filamentous Algae Percent		Moss Percent		Wood Cover - Unembedded		Wood Cover - Total		Round Rock Cover - Unembedded		Round Rock Cover - Total		Flat Rock Cover - Unembedded		Flat Rock Cover - Total
			Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.
fixed cond	(Intercept)		1.866***	0.386	5.699***	0.252	2.762***	0.320	5.417***	0.381	2.373***	0.355	2.580***	0.453	2.407***	0.523	5.441***	0.263	-0.527**	0.190	-0.204	0.202	-3.971***	0.690	-3.585***	0.664	-3.139***	0.653	-1.978***	0.462	-1.736***	0.473	-2.585***	0.562	-2.529***	0.638	-4.790***	0.712	-4.905***	0.792
	`Upstream catchment developed (%)`		0.002	0.006	-0.004	0.004	0.005	0.005	-0.001	0.006	0.007	0.005	0.017*	0.007	0.023**	0.008	-0.010*	0.004	0.009**	0.003	-0.003	0.003	0.035***	0.011	0.021*	0.009	0.018+	0.010	-0.015*	0.007	-0.018*	0.007	0.011	0.009	0.013	0.010	0.030**	0.011	0.033**	0.012
	`Upstream catchment agricultural (%)`		-0.006	0.006	-0.004	0.004	0.004	0.005	-0.006	0.006	0.006	0.005	0.012+	0.007	0.020*	0.008	-0.013**	0.004	0.011***	0.003	-0.004	0.003	0.034**	0.010	0.016+	0.009	0.019+	0.010	-0.017*	0.007	-0.020**	0.007	0.005	0.008	0.004	0.009	0.032**	0.011	0.034**	0.012
	log1p(`Reach Slope`)		-0.317	0.237	-0.579***	0.152	0.497*	0.195	-0.140	0.233	0.005	0.192	1.157***	0.281	1.086***	0.325	-0.011	0.162	0.126	0.113	0.145	0.117	0.607	0.390	0.852*	0.370	1.156**	0.405	0.436	0.289	0.500+	0.291	1.054**	0.345	1.327***	0.393	1.643***	0.441	1.930***	0.493
	Year_n		0.001	0.006	0.001	0.007	0.003	0.007	0.020**	0.008	0.011	0.013	-0.008	0.012	-0.003	0.010	-0.001	0.009	0.004	0.006	-0.004	0.007	-0.021	0.029	0.016	0.026	-0.008	0.023	0.031*	0.015	0.024	0.015	-0.025	0.017	-0.028+	0.017	0.024	0.018	0.019	0.019
	`UCA Class`Intermediate		0.397***	0.115	0.133+	0.076	0.310**	0.098	0.637***	0.121	0.031	0.128	0.334*	0.150	0.171	0.158	-0.051	0.095	-0.084	0.067	0.013	0.079	0.041	0.209	0.309	0.233	0.293	0.236	0.298	0.214	0.287	0.220	0.062	0.177	0.168	0.208	0.285	0.220	0.406+	0.239
	`UCA Class`Large		0.829***	0.141	0.267**	0.092	0.564***	0.119	1.344***	0.148	-0.007	0.140	0.556**	0.181	0.396*	0.190	-0.108	0.115	0.111	0.082	0.036	0.094	0.209	0.250	0.504+	0.259	0.775**	0.286	0.093	0.313	0.179	0.317	-0.363	0.222	-0.284	0.255	0.828**	0.263	1.003***	0.288
	`Habitat Stability Class`High		-0.414***	0.094	0.105+	0.062	-0.233**	0.080	-0.375***	0.101	-0.246*	0.105	-0.232+	0.126	-0.323*	0.129	0.214**	0.081	-0.091	0.058	-0.566***	0.068	-0.202	0.178	-0.293	0.188	-0.320	0.200	0.171	0.191	0.218	0.197	-0.288+	0.149	-0.408*	0.174	-0.367*	0.178	-0.395*	0.192
	`Stress Class`Medium		0.130	0.194	0.001	0.126	-0.157	0.162	0.166	0.198	0.471*	0.196	0.275	0.246	0.054	0.267	0.241	0.148	-0.229*	0.103	-0.147	0.113	-0.152	0.368	0.017	0.403	0.789*	0.366	0.037	0.264	-0.035	0.273	0.354	0.306	0.475	0.353	0.262	0.377	0.419	0.414
	`Stress Class`High		-0.114	0.309	0.165	0.200	-0.115	0.260	0.322	0.312	0.188	0.281	-0.209	0.373	-0.521	0.431	0.217	0.219	-0.226	0.155	-0.134	0.164	-0.937+	0.545	-0.267	0.545	0.197	0.557	0.459	0.429	0.419	0.433	-0.172	0.468	-0.131	0.532	-0.419	0.595	-0.348	0.666
	Year_n × `UCA Class`Intermediate		-0.001	0.007	0.001	0.008	-0.005	0.008	-0.008	0.009	0.017	0.013	0.017	0.014	0.002	0.012	0.013	0.010	0.002	0.007	0.005	0.008	0.040	0.030	0.044+	0.025	-0.009	0.024	-0.008	0.016	-0.006	0.016	0.035+	0.020	0.026	0.019	-0.012	0.022	-0.008	0.023
	Year_n × `UCA Class`Large		-0.003	0.008	0.006	0.010	-0.018+	0.010	-0.019	0.012	0.014	0.018	0.006	0.019	-0.008	0.016	0.009	0.014	-0.021*	0.010	0.007	0.011	0.031	0.035	0.019	0.035	-0.047	0.048	-0.011	0.028	-0.009	0.028	0.039	0.029	0.033	0.027	-0.029	0.031	-0.028	0.031
	Year_n × `Habitat Stability Class`High		0.000	0.004	-0.002	0.005	-0.008+	0.005	-0.006	0.005	0.007	0.008	0.007	0.008	-0.000	0.006	0.008	0.007	0.001	0.005	0.005	0.005	0.002	0.014	0.008	0.015	0.002	0.016	0.001	0.013	-0.001	0.013	-0.012	0.010	-0.008	0.010	-0.014	0.012	-0.009	0.012
	Year_n × `Stress Class`Medium		0.003	0.006	-0.001	0.007	0.006	0.007	-0.008	0.008	-0.013	0.013	0.009	0.013	0.001	0.011	-0.007	0.009	0.003	0.007	0.016*	0.007	0.047	0.028	0.041	0.026	-0.049*	0.022	-0.010	0.013	0.001	0.014	0.028	0.018	0.021	0.017	0.002	0.019	-0.001	0.019
	Year_n × `Stress Class`High		0.005	0.006	0.007	0.006	0.002	0.007	-0.005	0.007	-0.006	0.013	0.015	0.012	0.006	0.010	-0.012	0.009	-0.004	0.006	0.012+	0.007	0.063*	0.028	0.032	0.025	-0.031	0.022	-0.031*	0.015	-0.021	0.015	0.044*	0.018	0.029+	0.017	0.029	0.018	0.020	0.018
	Year_n × `UCA Class`Intermediate × `Stress Class`Medium		-0.002	0.008	0.008	0.009	-0.009	0.009	0.004	0.011	-0.016	0.014	-0.029+	0.016	-0.005	0.014	-0.009	0.011	0.003	0.008	-0.012	0.009	-0.047	0.031	-0.043	0.027	-0.001	0.026	0.002	0.019	0.008	0.018	-0.027	0.022	-0.024	0.021	-0.017	0.025	-0.026	0.026
	Year_n × `UCA Class`Large × `Stress Class`Medium		0.005	0.009	-0.004	0.010	0.009	0.011	0.013	0.013	-0.011	0.018	-0.013	0.019	0.010	0.017	-0.009	0.014	0.015	0.010	-0.012	0.010	0.010	0.035	0.002	0.035	0.075	0.049	-0.003	0.030	-0.007	0.030	-0.022	0.030	-0.016	0.028	0.014	0.032	0.018	0.032
	Year_n × `UCA Class`Intermediate × `Stress Class`High		-0.000	0.008	0.001	0.008	-0.011	0.009	0.005	0.010	-0.009	0.014	-0.001	0.015	0.012	0.014	-0.009	0.011	0.003	0.008	-0.000	0.008	-0.027	0.030	-0.025	0.027	0.014	0.026	0.001	0.018	0.003	0.018	-0.014	0.022	-0.006	0.021	0.007	0.024	0.007	0.026
	Year_n × `UCA Class`Large × `Stress Class`High		0.001	0.009	-0.012	0.010	0.010	0.011	0.001	0.013	0.002	0.018	-0.012	0.019	0.006	0.017	-0.006	0.014	0.015	0.010	0.001	0.011	-0.019	0.036	-0.012	0.036	0.033	0.048	-0.003	0.030	-0.009	0.030	-0.032	0.031	-0.018	0.029	-0.005	0.032	0.011	0.033
	Year_n × `UCA Class`Intermediate × `Habitat Stability Class`High		0.004	0.006	0.001	0.007	0.011	0.007	0.004	0.009	0.001	0.010	0.004	0.011	-0.001	0.009	-0.003	0.009	-0.003	0.007	0.003	0.007	0.017	0.020	-0.014	0.020	0.013	0.020	-0.006	0.017	-0.012	0.017	-0.002	0.014	0.003	0.014	0.026	0.017	0.023	0.018
	Year_n × `UCA Class`Large × `Habitat Stability Class`High		-0.001	0.006	0.011+	0.007	0.006	0.007	0.022**	0.008	0.004	0.009	-0.003	0.011	0.002	0.009	-0.001	0.009	0.009	0.007	0.004	0.007	-0.041*	0.020	-0.027	0.020	-0.043+	0.026	0.022	0.025	0.017	0.025	0.000	0.015	0.000	0.015	0.019	0.018	0.013	0.018
ran_pars cond	sd__(Intercept)	`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode	0.455		0.253		0.344		0.483		0.241		0.495		0.580		0.218		0.211		0.183		0.393		0.234		0.631		0.273		0.351		0.494		0.677		0.707		0.790
	sd__Year_n		0.009		0.008		0.005		0.016		0.001		0.014		0.014		0.006		0.012		0.008		0.011		0.017		0.010		0.002		0.010		0.001		0.005		0.019		0.019
	cor__(Intercept).Year_n		-0.359		-0.061		-0.139		-0.496		-1.000		-0.435		-0.344		-0.132		-0.546		-0.562		1.000		1.000		-1.000		1.000		-0.426		1.000		-1.000		-0.178		-0.084
fixed zi	(Intercept)										9.441***	2.686	13.517*	5.544	16.431*	6.781	-6.457	11.758			-4.716	5.008	3.939***	1.171	5.023***	1.132	-2.260	1.747	-5.122**	1.778	-5.938**	1.845	4.733	4.058	7.115	4.360	11.185*	4.454	12.325+	6.685
	`Upstream catchment developed (%)`										-0.102*	0.041	-0.177*	0.087	-0.139	0.103	-0.001	0.116			0.164**	0.052	-0.027	0.018	-0.025	0.016	0.056*	0.026	0.073**	0.028	0.066*	0.028	0.010	0.056	-0.038	0.064	-0.052	0.064	0.020	0.098
	`Upstream catchment agricultural (%)`										-0.069+	0.039	-0.156+	0.080	-0.211*	0.104	0.011	0.121			0.207***	0.054	-0.030+	0.017	-0.016	0.015	0.051*	0.026	0.060*	0.027	0.057*	0.028	0.007	0.056	-0.012	0.064	-0.093	0.064	-0.029	0.095
	log1p(`Reach Slope`)										-7.655***	1.672	-11.966**	3.943	-13.408**	4.853	-16.559*	7.037			-5.617**	1.832	-1.563*	0.677	-2.305***	0.628	-3.582***	1.004	0.847	1.030	1.771+	1.071	-9.593***	2.703	-13.082***	3.119	-10.899***	3.001	-17.946***	4.647
	Year_n										0.056	0.067	-0.056	0.141	-0.013	0.190	0.733+	0.418			-2.753	1.708	-0.009	0.043	-0.144***	0.043	0.079	0.054	-0.053	0.054	-0.052	0.055	0.079	0.087	0.172	0.109	-0.026	0.106	-0.089	0.147
	`UCA Class`Intermediate										-1.609+	0.875	-3.283	2.060	-4.156*	1.874	6.689	8.824			-1.605	1.025	-1.185**	0.430	-1.094*	0.455	-2.461**	0.757	1.609**	0.622	1.926**	0.671	-3.129**	1.144	-3.828*	1.526	-4.845**	1.602	-8.342**	2.895
	`UCA Class`Large										-4.535***	1.203	-5.605*	2.781	-9.688*	4.052	-7.298	678.340			1.017	1.431	-0.921+	0.524	-1.540**	0.554	-2.015*	0.875	1.731*	0.818	2.264*	0.881	-1.668	1.204	-3.662*	1.752	-2.951	1.841	-5.247+	3.120
	`Habitat Stability Class`High										1.563*	0.793	2.778	2.033	2.705	2.278	6.581	7.414			-1.008	3.153	0.052	0.371	0.135	0.399	1.038	0.653	-0.383	0.532	-0.046	0.574	0.599	0.809	-0.355	1.225	1.200	1.250	0.894	2.094
	`Stress Class`Medium										-0.805	1.366	-1.277	2.748	-1.295	3.481	-7.097	9.221			-9.026**	3.221	-0.424	0.638	-1.134+	0.675	-2.726*	1.088	1.132	0.946	1.358	1.018	-3.042*	1.543	-2.906	2.171	-3.478	2.227	-7.985*	3.531
	`Stress Class`High										1.980	2.144	1.864	4.462	-2.495	5.080	-1.653	7.585			-9.511**	3.651	-0.385	0.947	-0.955	0.942	-3.028*	1.478	0.183	1.475	0.788	1.543	-3.111	2.447	-1.104	3.299	-4.915	3.381	-12.608*	5.389
	Year_n × `UCA Class`Intermediate										-0.084	0.079	-0.021	0.151	0.066	0.198	-0.770+	0.421			-11.963	38.848	-0.067	0.049	-0.006	0.044	0.059	0.056	-0.035	0.062	-0.025	0.061	-0.047	0.101	-0.120	0.116	0.052	0.115	0.181	0.154
	Year_n × `UCA Class`Large										0.062	0.110	0.064	0.185	-13.776	6240.302	-2.867	1074.702			-13.963	50.033	-0.146*	0.072	0.007	0.059	0.115	0.082	0.035	0.086	0.044	0.085	-0.134	0.131	-0.162	0.139	-0.092	0.147	-0.280	0.221
	Year_n × `Habitat Stability Class`High										-0.008	0.054	0.023	0.121	0.092	0.188	-0.384	0.334			0.958	1.236	0.001	0.034	0.001	0.033	-0.034	0.039	-0.002	0.038	-0.020	0.037	0.001	0.071	-0.034	0.082	-0.103	0.078	-0.084	0.109
	Year_n × `Stress Class`Medium										-0.064	0.068	-0.141	0.160	-0.285	0.216	0.041	0.344			1.686	1.314	-0.098*	0.046	-0.012	0.043	0.043	0.055	-0.074	0.056	-0.046	0.057	-0.118	0.083	-0.175+	0.106	0.140	0.099	0.230+	0.135
	Year_n × `Stress Class`High										-0.067	0.067	0.023	0.099	-0.217	0.177	-0.248	0.264			1.733	1.317	-0.073+	0.043	0.022	0.042	0.067	0.054	-0.005	0.056	0.002	0.058	-0.148+	0.081	-0.154	0.097	0.113	0.096	0.271*	0.136
	Year_n × `UCA Class`Intermediate × `Stress Class`Medium										0.026	0.092	0.287	0.216	0.238	0.264	-2.224	201.469			8.768	34.298	0.069	0.059	0.036	0.048	-0.023	0.058	-0.028	0.071	-0.035	0.071	-0.008	0.134	0.119	0.143	-0.073	0.138	-0.138	0.165
	Year_n × `UCA Class`Large × `Stress Class`Medium										-0.025	0.116	-3.081	42.130	7.183	6387.932	2.577	1055.116			0.027	137.236	0.112	0.075	0.020	0.058	0.102	0.081	0.082	0.091	0.042	0.089	-0.001	0.139	0.111	0.145	-0.321	0.222	-0.329	0.339
	Year_n × `UCA Class`Intermediate × `Stress Class`High										-0.053	0.098	-6.787	173.940	0.073	0.267	-6.192	452.218			8.648	34.298	0.012	0.061	-0.030	0.053	-0.128*	0.061	-0.055	0.072	-0.040	0.070	-0.020	0.141	0.025	0.156	-0.029	0.132	-0.056	0.162
	Year_n × `UCA Class`Large × `Stress Class`High										-0.068	0.127	-6.057	117.107	-0.146	16296.960	1.539	1337.650			0.632	137.236	0.103	0.077	-0.021	0.065	-0.031	0.081	-0.072	0.091	-0.090	0.089	0.115	0.137	0.085	0.142	-0.054	0.165	0.280	0.217
	Year_n × `UCA Class`Intermediate × `Habitat Stability Class`High										0.110	0.078	-0.039	0.187	-0.044	0.250	-25.059	10492.544			3.268	22.681	0.063	0.050	0.037	0.041	0.068	0.045	0.023	0.054	-0.010	0.052	0.065	0.125	0.074	0.124	0.107	0.112	0.111	0.137
	Year_n × `UCA Class`Large × `Habitat Stability Class`High										0.028	0.078	3.147	42.130	7.022	1365.423	-1.482	759.840			13.258	127.843	0.008	0.053	-0.005	0.045	-0.094+	0.051	-0.017	0.061	-0.011	0.059	0.081	0.105	0.151	0.111	0.301	0.185	0.236	0.249
ran_pars zi	sd__(Intercept)	`UCA Class`:`Stress Class`:`Habitat Stability Class`:SiteCode									2.446		3.349		3.831		7.023				1.294		0.603		0.896		2.174		1.385		1.708		2.211		4.021		3.955		7.444
	sd__Year_n										0.054		0.022		0.155		0.189				0.005		0.017		0.022		0.066		0.010		0.009		0.044		0.066		0.091		0.204
	cor__(Intercept).Year_n										-0.337		0.989		0.663		-0.997				0.995		1.000		-1.000		-0.906		1.000		-1.000		1.000		-0.591		-0.548		-0.924
	Num.Obs.		997		997		997		988		933		924		930		995		1006		952		997		997		997		997		997		997		997		997		997
	RMSE		3.50		62.28		30.57		89.28		15.79		64.59		43.46		49.69		0.09		0.11		0.18		0.18		0.15		0.04		0.04		0.10		0.10		0.10		0.11
	Pseudo R2		0.834		0.772		0.573		0.931		0.313		0.761		0.684		0.299		0.478		0.546		0.418		0.397		0.266		0.172		0.200		0.575		0.658		0.725		0.760
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;

Model Residuals

Results of the DHARMa package which uses a simulation-based approach to create readily interpretable scaled (quantile) residuals for fitted generalized linear (mixed) models. Left panel is a qq-plot to detect overall deviations from the expected distribution, with added tests for correct distribution (KS test), dispersion and outliers. Note that outliers are values that are defined as values outside the simulation envelope, not in terms of a particular quantile. Right panel is a plot of the residuals against the predicted value. Simulation outliers (data points that are outside the range of simulated values) are highlighted as red stars.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Model Fit

These figures show 1:1 plots of observed and model fitted values. The plots are shown on transformed scales. Endpoints were model fitted values fall more closely to the 1:1 line better fit the raw data. A good fitting model means the endpoint align with the distributional family in both mean and variance, that means of monitoring locations are normally distributed (on the transformed scale), and the overdispersion and zero-inflation terms are able to describe the attributes of the raw data. A pseudo R2 value is shown as represented by the squared spearman correlation coefficient of observed vs predicted values.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Results - Spatial Figures

Overall

This figure shows the global means and associated error bars. The Overall mean reflects the patterns across all monitoring locations, whereas the UCA Class, Stress Class, and Habitat Stability Class each reflect the patterns across locations assigned to those levels. For more detailed breakdowns across combinations of these 3 location classes, see the UCA, Stress and Stability tab.

UCA, Stress and Stability

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Endpoints by Monitoring Locations

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Spatial Correlations

Heatmap of spearman correlation coefficients among habitat endpoints (including Aquatic Ecosystem Classification variables), with green and blue indicating positive and negative relationships respectively. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Results - Temporal Figures

Overall

This figure shows the global temporal trend coefficients and associated error bars. The Overall trend reflects the patterns across all monitoring locations, whereas the UCA Class, Stress Class, and Habitat Stability Class each reflect the patterns across locations assigned to those levels. For more detailed breakdowns across combinations of these 3 location classes, see the UCA, Stress and Stability tab.

UCA, Stress and Stability

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Endpoints by Monitoring Locations

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Trend Correlation

Heatmap and cluster diagram of spearman correlations among temporal trends in habitat endpoints. Trends were coded as +2, +1, 0, -1, and -2 corresponding with ++, +, 0, -, and -- from the Trends Table. Strong positive correlations among habitat endpoints indicate shared increasing or decreasing temporal trends across multiple sites, while strong negative correlations indicate opposing increasing and decreasing temporal trends across multiple sites. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Results - Temporal Tables

Trends by Stream Order, Stress, and Stability

This table shows concurrent patterns of temporal trends across all habitat endpoints by stream order and stress class. The ++ and -- indicate that the 90 percent confidence interval of the slope coefficient does not intersect 0, suggesting a significant increasing and decreasing patterns, respectively. The + and - indicate that the 50 percent confidence interval of the of the slope coefficient does not intersect 0, suggesting increasing and decreasing non-significant trends.

Overall Trends	Morphology				Substrate Sizes			Indices			Percent			Cover
	Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Global Trend	++	++	-	++	++	+	0	0	+	++	++	++	--	+	+	+	0	++	++
Habitat Stability Class
Low	+	++	-	++	++	0	0	-	+	++	++	++	--	+	+	++	0	++	++
High	+	++	-	++	++	+	0	+	+	++	++	++	--	+	+	0	-	++	++
Stress Class
Low	0	+	-	++	++	0	-	+	0	+	0	+	-	++	+	0	-	+	+
Medium	++	++	0	++	+	0	0	-	++	++	++	++	--	+	++	+	-	++	+
High	++	++	-	++	++	++	+	-	0	++	++	++	--	0	0	++	+	++	++
UCA Class
Small	+	+	0	++	+	0	0	0	+	++	+	++	--	++	++	-	--	++	++
Intermediate	+	++	-	++	++	++	+	+	++	++	++	++	--	+	+	++	0	++	++
Large	+	++	-	++	++	0	0	0	0	++	+	++	--	+	0	+	0	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

UCA Class / Stress Class	Morphology				Substrate Sizes			Indices			Percent			Cover
	Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stress
UCA Small	++	++	0	++	+	+	0	-	0	++	++	++	--	0	0	+	0	++	++
UCA Intermediate	++	++	-	++	++	++	++	-	+	++	++	++	-	0	0	++	++	++	++
UCA Large	+	+	-	+	++	0	0	-	0	++	++	++	--	0	0	+	+	+	+
Medium Stress
UCA Small	+	0	+	+	0	0	0	0	+	++	++	++	--	+	++	0	-	+	+
UCA Intermediate	+	++	-	+	0	-	-	0	++	+	++	++	--	+	+	0	-	0	-
UCA Large	++	++	0	++	+	-	0	-	+	++	++	++	--	+	+	+	+	+	+
Low Stress
UCA Small	0	0	0	++	+	0	0	0	+	0	-	+	0	++	+	-	--	+	+
UCA Intermediate	0	0	0	+	++	+	0	+	+	+	+	++	0	+	+	0	0	+	+
UCA Large	0	+	-	+	+	0	-	+	-	+	0	+	-	+	+	0	0	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

UCA Class / Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
	Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stability
UCA Small	+	0	0	++	+	+	0	0	+	++	+	++	--	++	++	-	--	++	++
UCA Intermediate	++	+	-	++	++	++	0	+	+	++	++	++	-	0	0	+	0	++	++
UCA Large	0	++	-	++	++	0	0	0	0	++	0	+	--	+	+	0	0	+	+
Low Stability
UCA Small	+	+	+	++	0	0	0	-	+	+	+	++	--	+	+	0	-	++	++
UCA Intermediate	0	+	-	++	+	+	+	0	+	+	++	++	--	+	+	++	0	+	+
UCA Large	+	+	-	0	++	0	0	0	-	+	++	++	--	0	0	+	+	+	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

Trends by Location

This table shows concurrent patterns of temporal trends across all habitat endpoints by monitoring locations. The ++ and -- indicate that the 90 percent confidence interval of the slope coefficient does not intersect 0, suggesting a significant increasing and decreasing patterns, respectively. The + and - indicate that the 50 percent confidence interval of the of the slope coefficient does not intersect 0, suggesting increasing and decreasing non-significant trends.

HUMBER
HUMBER	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
HU008WM	Small	High	+	+	0	+	++	+	++	0	-	+	++	++	--	0	0	+	-	++	++
HU005WM	Small	High	0	+	-	+	++	+	0	0	-	+	++	++	--	0	0	0	0	+	+
HU006WM	Intermediate	High	--	0	-	-	++	++	++	0	0	++	++	++	-	0	0	++	++	++	++
Medium Stess Class
HU039WM	Small	High	0	0	0	0	+	0	0	0	+	+	+	++	--	++	+	-	-	0	0
HU016WM	Small	Low	+	0	+	++	0	0	0	-	+	+	++	++	--	+	+	0	0	+	0
HU001WM	Intermediate	High	0	+	0	-	+	0	0	0	0	+	++	++	--	+	+	0	0	0	0
HU028WM	Intermediate	High	0	+	0	-	+	0	0	0	+	+	++	++	--	+	+	0	-	0	0
HU026WM	Intermediate	High	+	+	0	+	+	0	0	0	0	+	++	++	--	+	+	0	-	0	0
HU027WM	Intermediate	High	+	0	0	+	+	0	0	0	+	+	++	++	--	+	+	0	--	0	0
HU011WM	Intermediate	High	+	0	0	0	0	0	0	0	0	++	++	++	--	+	+	0	-	0	0
HU015WM	Intermediate	High	0	0	0	0	+	-	0	0	+	+	++	++	--	+	+	0	0	0	0
HU014WM	Intermediate	High	0	0	0	0	+	0	0	0	0	+	++	++	--	+	+	0	-	0	0
HU022WM	Large	High	+	+	0	+	+	0	0	0	+	++	+	++	--	+	+	0	0	0	0
HU024WM	Large	High	+	+	0	++	+	0	0	0	0	+	+	++	--	++	+	0	0	0	0
HU013WM	Large	High	0	+	0	+	+	0	0	0	+	+	++	++	--	+	+	0	0	+	+
HU010WM	Large	Low	+	0	0	0	0	0	0	-	0	0	++	++	--	0	0	++	0	0	0
HU003WM	Large	Low	0	0	0	0	0	-	0	-	-	+	++	++	--	0	0	++	+	+	+
HU007WM	Large	Low	0	0	0	0	0	+	-	-	--	0	++	++	--	0	0	++	+	+	+
HU012WM	Large	Low	0	0	0	--	0	0	0	0	-	0	++	++	-	0	0	++	+	+	+
HU009WM	Large	Low	++	0	0	0	0	-	0	-	+	+	++	++	-	0	0	+	++	0	0
Low Stess Class
HU032WM	Small	High	0	0	0	+	+	0	0	0	+	0	0	+	0	++	+	--	--	0	0
HU030WM	Small	High	0	0	0	+	+	0	0	+	0	0	0	0	0	++	+	--	--	0	0
HU035WM	Small	High	0	0	0	+	+	0	0	0	+	0	0	+	0	++	+	--	--	0	0
HU033WM	Intermediate	High	0	0	0	0	++	0	0	+	0	+	++	++	0	+	0	0	0	+	+
HU036WM	Intermediate	Low	0	0	0	0	++	0	0	+	+	0	+	++	-	++	+	+	0	0	0
HU031WM	Intermediate	Low	-	0	0	++	++	0	0	+	++	0	+	++	-	++	+	+	0	0	0
HU037WM	Intermediate	Low	+	+	0	+	++	0	0	+	+	0	+	++	-	+	+	+	0	0	0
HU034WM	Intermediate	Low	0	0	0	+	++	+	0	+	0	0	+	++	-	++	+	0	0	0	0
HU038WM	Intermediate	Low	0	0	0	+	++	0	--	+	-	0	0	++	-	++	+	+	0	0	0
HU029WM	Large	Low	-	+	--	0	+	0	-	0	-	0	0	+	-	+	0	0	0	0	0
HU002WM	Large	Low	0	+	-	-	+	0	0	+	--	0	0	+	-	+	0	0	0	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

CARRUTHERS
CARRUTHERS	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Medium Stess Class
CC003WM	Small	High	0	0	0	+	+	0	-	0	0	++	+	++	--	++	+	-	-	0	0
CC002WM	Intermediate	High	++	++	0	++	+	0	-	0	+	+	++	++	--	+	+	0	0	0	0
CC001WM	Intermediate	High	+	+	0	++	+	0	0	0	+	+	++	++	--	+	+	0	0	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

HIGHLAND
HIGHLAND	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
HL001WM	Small	High	+	+	-	--	++	+	0	0	0	+	++	++	--	0	0	0	0	++	++
HL008WM	Small	High	+	0	0	++	++	+	0	-	-	+	++	++	--	0	0	+	-	++	++
HL009WM	Small	High	-	0	-	-	+	0	+	0	+	+	++	++	--	0	0	+	-	++	++
HL006WM	Small	Low	0	0	+	0	0	0	0	-	+	+	++	++	--	0	0	++	0	++	++
HL007WM	Small	Low	+	+	0	0	0	0	0	-	+	+	++	+	-	0	0	+	+	++	+
HL005WM	Intermediate	High	++	0	-	++	++	++	0	0	++	++	++	++	-	0	0	++	+	++	++
HL003WM	Intermediate	High	++	++	-	+	++	++	+	0	0	+	++	++	-	0	0	++	+	++	++
HL011WM	Intermediate	High	+	0	-	0	++	++	+	0	0	++	++	++	-	0	0	++	++	++	++
HL010WM	Intermediate	Low	++	++	-	++	+	+	+	0	0	+	++	++	-	0	0	++	++	++	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

DUFFINS
DUFFINS	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
DF001WM	Small	Low	-	+	0	0	0	0	0	-	-	0	++	+	--	0	0	++	0	++	++
DF002WM	Small	Low	+	0	+	+	0	0	0	-	0	+	++	++	--	0	0	++	0	++	++
Medium Stess Class
DF018WM	Small	High	0	0	0	0	+	0	0	0	+	++	+	++	--	++	+	-	-	0	0
DF005WM	Small	High	0	0	0	+	+	0	0	0	+	+	++	++	--	++	+	-	-	+	0
DF014WM	Small	High	0	0	0	-	0	+	+	0	0	+	++	++	--	++	+	-	--	+	+
DF016WM	Small	High	0	0	0	0	0	+	0	0	+	+	++	++	--	++	+	-	--	0	+
DF008WM	Small	Low	+	0	+	++	0	0	-	0	0	+	0	++	--	+	+	0	0	+	0
GTUR	Small	Low	0	0	+	0	0	0	0	-	0	+	+	++	--	+	+	0	0	+	0
DF007WM	Small	Low	+	0	+	+	0	0	0	-	0	+	+	++	--	+	+	0	-	+	+
DF020WMB	Intermediate	High	0	0	0	0	+	0	+	0	++	+	++	++	--	+	+	0	-	+	0
DF013WM	Intermediate	Low	0	0	-	0	0	0	+	0	++	+	+	++	--	+	+	+	0	0	0
DF011WM	Large	High	0	+	0	+	+	0	0	0	++	++	++	++	--	++	+	0	0	+	+
DF009WM	Large	High	++	++	0	++	+	0	-	0	+	+	++	++	--	++	+	0	0	+	0
DF010WM	Large	Low	0	0	0	0	0	-	--	-	+	+	++	++	-	0	0	++	+	0	0
DF003WM	Large	Low	0	0	0	++	0	0	+	-	0	+	++	++	-	0	0	++	+	+	0
Low Stess Class
DF021WMB	Small	Low	0	0	0	+	+	0	0	0	+	0	0	+	0	+	+	-	-	+	+
DF004WM	Small	Low	0	0	0	+	+	0	0	0	0	0	-	+	0	++	+	-	--	+	+
DF019WM	Intermediate	Low	0	0	0	+	++	0	0	+	+	0	+	++	-	++	+	+	0	0	0
DF017WM	Intermediate	Low	-	0	0	0	++	0	+	+	0	0	+	++	0	++	+	0	0	0	0
DF015WM	Intermediate	Low	0	0	0	+	++	+	0	+	+	0	+	++	-	+	+	+	0	+	+
DF012WM	Intermediate	Low	0	0	0	0	++	+	+	+	+	0	++	++	-	++	+	+	-	+	+
DF006WM	Large	Low	0	0	-	+	+	0	0	0	0	0	0	+	-	+	0	0	0	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

ETOBICOKE
ETOBICOKE	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
EC010WM	Small	High	+	+	0	++	++	+	0	0	0	+	++	++	--	0	0	+	-	++	++
EC006WM	Small	High	+	+	0	++	++	0	0	0	0	+	++	++	--	0	0	+	-	+	0
EC003WM	Small	High	+	+	0	++	++	0	--	0	-	+	++	++	--	0	0	0	-	+	0
EC009WM	Intermediate	High	-	+	-	-	++	++	++	0	0	+	++	++	-	0	0	++	+	++	++
EC001WM	Large	High	+	++	-	0	++	0	0	0	0	++	+	++	--	0	0	+	0	+	+
EC005WM	Large	Low	+	0	0	0	++	0	0	0	0	+	++	++	--	0	0	++	+	0	0
EC004WM	Large	Low	++	-	0	0	++	0	0	0	-	++	++	++	--	0	0	++	++	+	+
EC002WM	Large	Low	--	+	-	--	++	0	-	-	0	++	++	++	--	0	0	++	+	+	+
Medium Stess Class
Mayfield3	Intermediate	High	+	+	0	+	+	0	0	0	0	+	++	++	--	+	+	0	-	0	0
EC012WM	Intermediate	High	0	0	0	0	+	0	-	0	+	+	++	++	--	+	+	0	-	0	0
Mayfield2	Intermediate	High	0	++	0	0	0	0	0	0	0	+	++	++	--	+	+	0	--	0	0
EC007WM	Large	High	+	+	0	+	+	--	-	0	+	+	++	++	--	++	+	0	+	0	0
EC008WM	Large	Low	+	0	0	+	0	-	0	0	+	+	++	++	-	0	0	++	+	+	0
EC011WM	Large	Low	+	0	0	0	0	0	0	-	0	+	++	++	-	0	0	++	+	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

FRENCHMANS BAY
FRENCHMANS BAY	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
FB001WM	Small	High	+	0	0	++	++	+	0	0	0	+	++	++	--	0	0	0	0	++	+
FB002WM	Small	High	+	+	0	+	++	++	+	0	0	+	++	++	--	0	0	+	-	++	++
FB004WM	Small	High	+	++	0	+	++	+	0	0	0	+	++	++	--	0	0	+	0	++	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

DON
DON	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
DN011WM	Small	High	+	+	0	+	++	0	--	0	0	+	++	++	--	0	0	+	-	++	++
DN009WMB	Small	High	0	0	0	0	++	+	+	-	0	+	++	++	--	0	0	+	-	++	++
DN007WM	Small	High	0	0	0	0	+	0	0	-	0	+	++	++	--	0	0	+	-	++	++
DN014WM	Small	High	0	++	-	0	++	++	0	-	0	+	++	++	--	0	0	+	-	++	+
DN020WM	Small	High	0	0	-	0	++	+	0	0	0	++	++	++	--	0	0	0	0	++	++
DN018WM	Small	High	+	0	0	0	++	+	0	0	0	+	++	++	--	0	0	0	0	++	+
DN010WMB	Small	Low	+	+	+	+	0	0	0	-	0	+	++	++	--	0	0	++	0	++	+
DN006WMB	Small	Low	0	++	0	0	0	0	0	-	0	+	++	+	--	0	0	++	0	++	++
DN013WM	Small	Low	+	0	+	++	0	+	++	-	-	+	++	++	--	0	0	++	0	++	++
DN003WM	Intermediate	High	+	++	-	+	++	++	++	0	0	++	++	++	-	0	0	++	+	++	++
DN008WM	Intermediate	Low	0	0	-	0	+	+	+	-	+	+	++	++	-	0	0	++	++	++	+
DN002WM	Intermediate	Low	0	+	-	++	+	0	+	0	+	+	++	++	--	0	0	++	+	++	++
DN016WM	Intermediate	Low	+	++	--	+	+	+	0	0	0	+	++	++	--	0	0	++	+	+	+
DN001WM	Large	Low	++	0	0	+	++	0	0	-	-	0	++	+	--	0	0	++	++	+	+
DN005WM	Large	Low	0	0	0	0	++	0	0	-	0	++	++	++	--	0	0	++	+	+	+
DN012WM	Large	Low	+	0	0	+	++	0	0	-	0	+	++	++	--	0	0	++	+	0	0
DN015WM	Large	Low	0	0	0	0	++	0	+	-	0	+	++	++	--	0	0	++	+	+	+
Medium Stess Class
DN019WM	Small	High	0	0	0	+	+	+	+	0	+	++	+	++	--	++	+	-	-	0	0
DN023WM	Small	Low	0	0	+	0	0	0	0	-	+	+	++	++	--	+	+	0	-	+	+
DN021WM	Intermediate	Low	+	+	0	+	0	-	0	0	+	0	+	++	--	+	+	+	0	0	0
Low Stess Class
DN022WM	Small	High	0	0	0	0	+	0	0	+	-	0	-	+	0	++	+	--	--	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

ROUGE
ROUGE	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
RG004WM	Small	High	0	0	0	0	++	+	0	0	+	+	++	++	--	0	0	+	-	++	+
RG005WM	Small	High	0	+	0	0	++	0	0	0	+	+	++	++	--	0	0	+	0	++	+
RG017WM	Small	High	+	+	0	++	++	0	0	0	0	+	++	++	--	0	0	+	-	+	+
RG024WM	Small	High	+	0	0	+	++	+	+	-	0	+	++	++	--	0	0	+	-	++	+
RG003WM	Small	Low	+	+	+	++	0	0	0	-	0	+	++	++	--	0	0	++	0	++	+
Medium Stess Class
RNUP-3	Small	High	0	0	0	0	+	0	0	0	0	++	++	++	--	++	+	-	-	0	0
RNUP-1	Small	High	0	-	0	0	+	0	0	0	+	++	+	++	--	++	+	-	-	0	0
RG022WM	Small	High	+	0	0	+	0	+	0	0	0	++	++	++	--	++	+	0	--	+	+
RG023WM	Small	High	+	0	0	+	0	0	0	0	+	++	++	++	--	++	+	0	--	+	+
RG010WM	Small	Low	0	0	+	0	0	0	0	0	+	+	++	++	--	+	+	0	-	++	+
RNUP-2	Intermediate	High	0	0	0	0	+	0	0	0	+	+	++	++	--	+	+	0	-	0	0
RG011WM	Intermediate	High	+	+	0	+	+	0	-	0	0	+	++	++	--	+	+	0	-	0	0
RG018WM	Intermediate	High	+	+	0	+	0	0	0	0	+	++	++	++	--	+	+	0	--	+	0
RG016WM	Intermediate	High	0	+	0	0	+	-	-	0	+	+	++	++	--	+	+	0	-	+	0
RG019WM	Intermediate	High	+	++	0	+	+	0	0	0	+	+	++	++	--	+	+	0	0	0	0
RG020WM	Intermediate	High	+	++	0	+	+	0	0	0	0	++	++	++	--	+	+	0	0	0	0
RG021WM	Intermediate	High	0	0	0	0	+	0	0	0	+	++	++	++	--	+	+	0	-	+	0
RG015WM	Intermediate	Low	0	+	0	+	0	-	0	0	+	0	+	++	--	+	+	+	0	0	-
RG007WM	Large	High	++	+	0	++	+	0	-	0	++	+	++	++	--	++	+	0	0	0	0
RG009WM	Large	High	0	++	0	++	+	+	++	0	+	++	+	++	--	+	+	0	0	+	+
RG001WM	Large	Low	+	0	0	+	0	0	++	0	++	+	++	++	-	0	0	++	+	0	0
RG006WM	Large	Low	+	0	0	0	0	0	0	-	0	0	++	++	-	0	0	++	+	0	0
RG002WM	Large	Low	+	+	0	0	0	-	+	-	0	+	++	++	-	0	0	++	0	+	0
Low Stess Class
RG026WM	Small	High	-	0	-	++	+	0	0	0	0	0	-	+	0	++	+	--	--	0	0
RG027WM	Small	High	0	0	0	+	+	0	0	+	+	0	0	+	0	++	+	--	--	0	0
RG025WM	Small	High	0	0	0	0	+	0	-	0	0	0	0	+	0	++	+	--	--	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

MIMICO
MIMICO	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
MM005WM	Small	High	+	+	0	+	++	++	+	0	0	++	++	++	--	0	0	+	-	++	++
MM004WM	Intermediate	High	+	++	-	+	++	++	+	0	0	++	++	++	-	0	0	++	++	++	++
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

PETTICOAT
PETTICOAT	UCA Class	Habitat Stability Class	Morphology				Substrate Sizes			Indices			Percent			Cover
			Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Medium Stess Class
PT004WM	Small	High	+	0	0	++	0	+	-	0	+	++	++	++	--	++	+	-	--	+	0
PT002WM	Intermediate	High	+	+	0	+	+	0	0	0	+	+	++	++	--	+	+	0	-	+	0
PT001WM	Intermediate	High	0	+	0	0	+	0	--	0	0	++	++	++	--	+	+	0	--	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

Results - Habitat - Biota Relationships

Correlations

Heatmap of spearman correlation coefficients of habitat and fish community endpoints, with green and blue indicating positive and negative relationships respectively. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Scatterplots

Scatterplots of habitat and fish community endpoints where absolute spearman correlation coefficients are greater than 0.3.

Site Specific Evaluations

By Les Stanfield

A subset of sites were identified which field crews had identified visually as becoming “more degraded over time”. The intention was to provide a quality control check on the broader analysis. That is, are the broad patterns in the analysis in agreement with field impressions. Such intuitive interpretation provides a form of quality control for the conclusions.

A Crew leader (Mark Szonda) and Jeff Vandenburg (program supervisor) developed the following list of sites deemed to be undergoing “significant morphological change over the observed time of sampling.

The selected sites were: * CC002WM * DN008WM * HU021WM (excluded due to lack of data) * HU032WM * DF005WM * EC002WM * HU024WM * PT001WM * DF008WM * EC012WM * HU029WM * RG003WM

Interpretations from the trends by location table were summarized to provide a synopsis of directions for each habitat endpoint. Each endpoint was classified as to whether over the 25 year period there was no change (0); non significant increasing (+1) or decreasing (-1) change; or significant increasing (+2) or decreasing (-2) change. Habitat variables were ranked according to the mean change over time for all the sites for that category. For example, a value close to 2 indicates that nearly all the sites had significant changes in that variable. The results of this analysis is summarized below:

TO UPDATE

This analysis supports the broader conclusions that the most consistent measure of change over time in habitat is reflected in an increase in the amount of algae, both non and filamentous with a concurrent decrease in the amount of moss at a site. These conclusions are strong and consistent with ecological theory. Degrading streams have less stable substrates to which moss would attach. Algae is a much faster and seasonal plant and would be more adapted to grow on less stable substrates. Additionally, degrading streams are often found to receive increased nutrients which would also encourage algal growth.

The amount of wood at a site also tended to increase. This tendency is more complicated to interpret. Wood would be hypothesized to increase under two scenarios. Stable streams with well forested riparian zones tend to have large volumes of wood. Wood tends to decrease when riparian corridors are converted from mature forests to meadows or scrublands. Finally, under a scenario of decreasing stability, an expectation would be that wood would be increasingly mobilized into the channel through erosion and or migration from upstream. This interpretation assumes greater stream instability that is not consistently supported by the various geomorphic metrics.

None of the geomorphic metrics consistently indicate a degradation in channel stability and in some cases the data are confounding. For example, in contrast to expectations, bank stability tended to increase at most of the sites over time. A cursory look at the 2025 site features data for these sites determined that at only one site (EC012WM) was there evidence of channel hardening, negating this activity as an explanation of these findings. There were also no consistent changes in the width or depth, substrate or sediment sorting. At many of the site’s crews recorded evidence of sediment loading that could be a factor in their classification of these sites as being degrading over time. Perhaps there is a more nuanced and multivariate relationship between these metrics and channel stability, but there is no simplistic correlation between what biologists interpret as channel degradation and the OSAP derived metrics of channel stability. We will continue to seek an independent measure of channel degradation and also to reevaluate the existing metrics of channel stability in attempts to improve our ability to measure this important stream metric. In the meantime, we have been able to conclude that moss really can’t grow on a rolling stone, and perhaps this insight could direct future analysis and interpretations.

Results and Conclusions

Model Fit

Overall models fit observed data well. Generally, the best fits were observed for Morphology endpoints, and weakest for some Substrate Sizes, Percent, and Cover endpoints. Most modeled endpoints contained some statistical deviations from expectations based on the selected distributional family, but visual inspection of residual plots showed good matches between observed and expected values to their theoretical quantiles.

Spatial Patterns

In FWIS, the habitat stability score is a calculation based on habitat homogeneity, width:depth ratios, sediment sorting index, and sediment transport potential. The first three of these are included as separate endpoints in these analysis, but each site was assigned only a single stability class (high or low). Therefore, the spatial patterns attributed to habitat stability in habitat homogeneity, width:depth ratios, and sediment sorting index are artifacts of these endpoints being used to calculate stability. The habitat stability predictor was included primarily to interpret temporal changes, but can provide some additional insights into spatial patterns not directly associated with the calculation of stability (i.e., percent and cover endpoints).

Spatial Correlations

At a high level, there are groups of endpoints with positive and negative correlations. Flat and round rock cover, substrate sizes, and percent algae and moss endpoints all showed some degree of positive correlations. Most of these endpoints were negatively correlated with wood cover and sediment sorting index. These endpoints were also positively correlated with reach slope, upstream catchment developed (%) and upstream catchment stress index, but uncorrelated with upstream agriculture. Stream morphology endpoints were also positively correlated with this cluster, but only weakly so.

Morphology

Overall there were clear and expected patterns in morphology across the Upstream Catchment Area (UCA) classes. For stream width, there was a non-significant pattern of high stability streams being narrower, especially among smaller streams. Conversely, for stream depth, there was a non-significant pattern of low stability streams being shallower, but only among bigger streams. These patterns were consistent regardless of stress class. Together these patterns showed higher usable areas among low stability streams, which were significant among high and medium stress small sites.

Substrate Sizes

Substrate sizes varied predictably with UCA, as smaller streams had smaller substrates. While substrate sizes are expected to vary by landscape stress and habitat stability, the models were unable to consistently attribute spatial differences in these endpoints to either of these predictors. This suggests that additional modifying factors not accounted for in these models may be driving spatial patterns in substrates. Many of the monitoring locations with the smallest substrate sizes were found in Carruthers Creek watershed, but also some among the smaller streams in the upper Don and Humber watersheds.

Indices

Overall mean bank stability was lowest in high stress streams, and highest in low stress streams. Within the low stress streams, there is a pattern of decreasing bank stability among larger streams than smaller streams. Otherwise, there are few strong spatial patterns associated with UCA or stress classes among sediment sorting index or habitat homogeneity..

Percent

Among the percent endpoints, there is a greater amount of variation among larger streams, than small and intermediate streams, and again greater variation among low stress streams than medium and high stress streams. This variation appears to be driven by habitat stability, with variation in stability seeming to impact large streams more than small and intermediate streams. One potential explanation for this pattern is that there is a threshold response to growth of algae and moss in relation to some driving factor such as stream flow, where a minimal threshold of stream flow is required for the effects of habitat stability to be observed on moss and algae formation. Filamentous and non-filamentous algae tended to be high in sites closer to lake Ontario and increasing further from the lake. Carruthers Creek and Frenchmans Bay watersheds locations were an exception and had high percent algae. Percent moss did not show this spatial pattern. High algae percentages may be related to elevated nutrient concentrations among more agrigultural streams, especially in the upper Humber river watershed.

Cover

There is little difference in the spatial patterns between the unembedded and total fractions of wood, round rock and flat rock cover endpoints. There are also few consistent patterns among the cover endpoints in relation to overall UCA or stress. Wood and flat rock cover especially didn’t show any significant patterns at higher levels. Round rock cover did show significant patterns where low stability sites had significantly higher round rock cover than high stability sites. Wood cover tended to be lowest closer to lake Ontario, and higher further from the lake. This may be reflective of reduced forest cover in more densely urbanized areas closer to the lake. Round and flat rock covers both tend to be higher closer to lake Ontario, suggesting greater proportions of these cover types in more developed areas. These spatial patterns are confirmed in trend analysis which showed positive correlations with Upstream Catchment Stress Index and Upstream Catchment Developed (%).

Temporal Patterns

Temporal Correlations

Temporal trends in flat and round rock cover, substrate sizes and stream morphology were generally positively correlated, indicating share increasing or decreasing trends over time. Wood cover, width:depth ratio, bank stability, sediment sorting index and moss percent were mostly negatively correlated with the previous endpoints, indicating contrasting trends (i.e., one endpoint is increasing while the other is decreasing over time, or vice versa). Wood cover, width:depth ratio, bank stability, sediment sorting index and moss percent were generally positively correlated with each other.

Morphology

Stream width, depth and usable area are showing overall patterns of increasing over time throughout the TRCA. The rate of increasing width and depth was greatest in high and medium stress streams. Habitat stability had little overall effect on temporal changes. Width to depth ratios were on average decreasing, but showed considerable and inconsistent variability between stream sizes, stress classes and habitat stability.

Substrate Sizes

The D50 endpoint showed consistently increasing trends over time, especially at large and intermediate streams. There are also patterns of increasing D50 at high and low stress classes, but not medium stress. High stability streams had slightly (but not significantly) greater rates of change than low stability streams. Variation in temporal trends among D85 and D50Max were inconsistent across UCA, stress, and stability classes. This suggests there are is little evidence temporal trends in these endpoints are related to broader landscape patterns.

Indices

Habitat homogeneity showed consistent increasing trends over time. The rate of increase was highest in high stress streams, followed by medium stress. Low stress streams showed no evidence of temporal changes in homogeneity. There were consistent, but not statistically significant differences in the rate of change between high and low stability streams, with high stability streams changing more rapidly. Bank stability and sediment sorting index did not show consistent or significant temporal changes across UCA, stress or habitat stability classes.

Percent

Filamentous and non-filamentous algae percents both increased over time at high and medium stress sites, and showed no consistent trends at low stress sites. The trends did not vary consistently by upstream catchment area, but among large, high-stress streams, sites with high habitat stability tended to show a higher rate of change than low stability sites. Rates of change were consistently higher for non-filamentous algae than for filamentous algae. Moss percent had an inverse pattern from algae, showing significant decreases over time at high and medium stress sites, and slightly faster rates of decline among large high stability sites than small high stability sites.

Cover

There is little difference in the temporal patterns between the unembedded and total fractions of wood, round rock and flat rock cover endpoints. There were slight overall increasing trends in wood cover, but only among small streams. Round rock cover showed slight increasing trends, but only among high stress streams, and slight decreases over time at small, low stress streams. Flat rock cover showed similar increasing trends across all sizes of high stress streams. Overall cover endpoints showed inconsistent patterns across size, stress and habitat stability classes.

Relationships Between Habitat and Fish Communities

There are no strong relationships (spearman rho > 0.6) between habitat and fish community endpoints. Some weak, but noteworthy relationships do exist however. Specifically, lithophilic spawners and open substratum spawners were positively correlated with a cluster of habitat endpoints including round and flat rock cover, substrate sizes, moss percent and habitat homogeneity. Sensitive darters were positively correlated with width, usable area, and width:depth ration, but not depth alone. Interestingly total community biomass, creek chub and community density were negatively correlated with these same morphology endpoints. Wood cover was also positively correlated with brook trout, creek chub, and coldwater stenotherm fish endpoints. Fish community richness, tolerant taxa richness and SATI biomass (a thermally weighted biomass measure) were positively correlated with width, depth and usable area.

Conclusions

Overall, morphology and percent endpoints showed strongest and most consistent spatial and temporal patterns to broad landscape level drivers. Morphology also showed some of the strongest correlations with fish community endpoints. Cover endpoints showed less consistent spatial and temporal patterns, but did show some relationships with fish community endpoints.

The stress classes used in the current models are perhaps suboptimal for describing spatial patterns in stream habitat data. Though correlation between most habitat endpoints were stronger with the Upstream Catchment Stress Index than with Upstream Catchment Developed (%), perhaps different thresholds should be considered to better understand spatial patterns among in stream habitat endpoints.